Setting the Working Directory and Global Code Chunk Options

Using include = FALSE means that this chunk will run, but the code and the results will not appear in the knitted document. We have set message and warning to FALSE for every chunk so the knitted document is not bombarded with messages and warnings.

Packages to Load

There is a brief explanation of each package’s purpose left as a comment following the library call.

Import Data

The document “webs” is a text file in the data folder of the repository/directory.

Data Wrangling

Here, we rename some levels of factors Location and Land and put them in the order we want the categories to appear in a graph. We also create another data set “sites” that restricts “webs” to just a single data point for each site. We use “sites” to compare search distance, the number of webs, and the number of spiders. We will use “webs” for web height. The “web_dist” data set is a subset of “webs” that looks at the distance from the focal web to the nearest and second nearest neighbors.

At one time, we tried scaling and centering the data in difference ways. I want to make mention of it here because it might be useful information for me to look back on in the future.

It’s important to remember to back transform when making figures. We didn’t scale and center because log-transformations on right-skewed variables were enough to remove convergence errors.

Transformation of Variables

Many data are right-skewed. We check the distribution of the data and check distributions after log-transformation for right-skewed data not containing zeros. Echo = FALSE shows the results, but does not show the code used.

## 
##  Shapiro-Wilk normality test
## 
## data:  sites$tree100m
## W = 0.79926, p-value = 0.0004871
## 
##  Shapiro-Wilk normality test
## 
## data:  sites$tree100m_frac
## W = 0.79926, p-value = 0.0004871

## 
##  Shapiro-Wilk normality test
## 
## data:  sites$imperv100m
## W = 0.74128, p-value = 6.883e-05

## 
##  Shapiro-Wilk normality test
## 
## data:  sites$TotalSub
## W = 0.85215, p-value = 0.003707
## 
##  Shapiro-Wilk normality test
## 
## data:  log(sites$TotalSub)
## W = 0.96252, p-value = 0.5417

## 
##  Shapiro-Wilk normality test
## 
## data:  sites$Traffic_Dist_Road
## W = 0.70845, p-value = 2.511e-05
## 
##  Shapiro-Wilk normality test
## 
## data:  log(sites$Traffic_Dist_Road)
## W = 0.93619, p-value = 0.1651

## 
##  Shapiro-Wilk normality test
## 
## data:  sites$Traffic_Dist_Highway
## W = 0.59575, p-value = 1.195e-06
## 
##  Shapiro-Wilk normality test
## 
## data:  log(sites$Traffic_Dist_Highway)
## W = 0.9091, p-value = 0.04539

## 
##  Shapiro-Wilk normality test
## 
## data:  sites$spec_rad
## W = 0.85174, p-value = 0.003645
## 
##  Shapiro-Wilk normality test
## 
## data:  log(sites$spec_rad)
## W = 0.84964, p-value = 0.003345

## 
##  Shapiro-Wilk normality test
## 
## data:  sites$light_rad
## W = 0.84052, p-value = 0.002319
## 
##  Shapiro-Wilk normality test
## 
## data:  log(sites$light_rad)
## W = 0.81646, p-value = 0.0009161

## 
##  Shapiro-Wilk normality test
## 
## data:  sites$patch_area_mm
## W = 0.70203, p-value = 2.077e-05
## 
##  Shapiro-Wilk normality test
## 
## data:  log(sites$patch_area_mm)
## W = 0.89878, p-value = 0.0281

## 
##  Shapiro-Wilk normality test
## 
## data:  sites$road_length_m
## W = 0.87899, p-value = 0.01156

## 
##  Shapiro-Wilk normality test
## 
## data:  sites$trail_length_m
## W = 0.7756, p-value = 0.0002127

We conclude to use log-transformations for TotalSub, Traffic_Dist_Road, Traffic_Dist_Highway, and patch_area_km.

Our analyses include looking at variation within each location. As such, this chunk subsets the “sites” and “webs” data by Location.

Variable Reduction Based on Correlations

Since many of these variables are likely highly correlated, we want to reduce the variables that we include in our model by removing variables that are highly correlated.

corr <- sites %>% 
  dplyr::select(tree100m, imperv100m, TotalSub, 
         Traffic_Dist_Road, Traffic_Dist_Highway, 
         spec_rad, light_rad, patch_area_mm, road_length_m)
# We exclude trail length from global correlations because trail length can only collected in Wilderness Park, and not for UNL City Campus

# When running findCorrelations, we get a vector of variables to remove to reduce pairwise correlations. 

findCorrelation(cor(corr, method = "spearman"), cutoff = .6, verbose = TRUE, names = TRUE) 
## Compare row 6  and column  2 with corr  0.73 
##   Means:  0.693 vs 0.594 so flagging column 6 
## Compare row 2  and column  8 with corr  0.77 
##   Means:  0.68 vs 0.568 so flagging column 2 
## Compare row 8  and column  9 with corr  0.738 
##   Means:  0.633 vs 0.53 so flagging column 8 
## Compare row 9  and column  7 with corr  0.706 
##   Means:  0.613 vs 0.487 so flagging column 9 
## Compare row 7  and column  1 with corr  0.623 
##   Means:  0.601 vs 0.44 so flagging column 7 
## Compare row 1  and column  3 with corr  0.662 
##   Means:  0.494 vs 0.336 so flagging column 1 
## All correlations <= 0.6
## [1] "spec_rad"      "imperv100m"    "patch_area_mm" "road_length_m"
## [5] "light_rad"     "tree100m"
# This suggests that we remove tree100m, imperv100, spec_rad, light_rad, patch_area_km, road_length_m
# That leaves us with TotalSub, Traffic_Dist_Road, and Traffic_Dist_Highway

chart.Correlation(corr, histogram = TRUE, method = "spearman") # kendall, spearman

corr <- sites %>% 
  dplyr::select(TotalSub, Traffic_Dist_Road, Traffic_Dist_Highway)
chart.Correlation(corr, histogram = TRUE, method = "spearman") # kendall, spearman

# Let's also test the variables after transformation

corr_transformed <- sites %>% 
  dplyr::select(tree100m, imperv100m, log_TotalSub, 
         log_tdr, log_tdh, 
         spec_rad, light_rad, log_patch, road_length_m)

findCorrelation(cor(corr_transformed, method = "spearman"), cutoff = .6, verbose = TRUE, names = TRUE) 
## Compare row 6  and column  2 with corr  0.73 
##   Means:  0.693 vs 0.594 so flagging column 6 
## Compare row 2  and column  8 with corr  0.77 
##   Means:  0.68 vs 0.568 so flagging column 2 
## Compare row 8  and column  9 with corr  0.738 
##   Means:  0.633 vs 0.53 so flagging column 8 
## Compare row 9  and column  7 with corr  0.706 
##   Means:  0.613 vs 0.487 so flagging column 9 
## Compare row 7  and column  1 with corr  0.623 
##   Means:  0.601 vs 0.44 so flagging column 7 
## Compare row 1  and column  3 with corr  0.662 
##   Means:  0.494 vs 0.336 so flagging column 1 
## All correlations <= 0.6
## [1] "spec_rad"      "imperv100m"    "log_patch"     "road_length_m"
## [5] "light_rad"     "tree100m"
# This suggests to remove tree100m, imperv100m, spec_rad, light_rad, log_patch, and road_length_m.

# This leaves us with log_TotalSub, log_tdr, and log_tdh, the same results as the non-transformed. 

chart.Correlation(corr_transformed, histogram = TRUE, method = "spearman") # kendall, spearman

corr_transformed <- sites %>% 
  dplyr::select(log_TotalSub, log_tdr, log_tdh)
chart.Correlation(corr_transformed, histogram = TRUE, method = "spearman") # kendall, spearman

corr_forest <- sites %>% 
  filter(Location == "Urban Forest") %>% 
  dplyr::select(tree100m, imperv100m, TotalSub, 
         Traffic_Dist_Road, Traffic_Dist_Highway, 
         spec_rad, light_rad, patch_area_km, road_length_m, trail_length_m)
# Notice trail length is included 

findCorrelation(cor(corr_forest, method = "spearman"), cutoff = .6, verbose = TRUE, names = TRUE) 
## Compare row 10  and column  2 with corr  0.696 
##   Means:  0.443 vs 0.315 so flagging column 10 
## Compare row 2  and column  9 with corr  0.899 
##   Means:  0.361 vs 0.286 so flagging column 2 
## Compare row 9  and column  1 with corr  0.679 
##   Means:  0.279 vs 0.272 so flagging column 9 
## Compare row 6  and column  7 with corr  0.685 
##   Means:  0.371 vs 0.258 so flagging column 6 
## Compare row 7  and column  3 with corr  0.911 
##   Means:  0.296 vs 0.219 so flagging column 7 
## All correlations <= 0.6
## [1] "trail_length_m" "imperv100m"     "road_length_m"  "spec_rad"      
## [5] "light_rad"
# This suggests removing imperv100m, spec_rad, light_rad, road_length_m, and trail_length_m
# We will keep tree100m, TotalSub, Traffic_Dist_Road, Traffic_Dist_Highway, and patch_area_km
chart.Correlation(corr_forest, histogram = TRUE, method = "spearman")

corr_forest <- sites %>% 
  filter(Location == "Urban Forest") %>% 
  dplyr::select(tree100m, TotalSub, Traffic_Dist_Road, patch_area_km, Traffic_Dist_Highway)
chart.Correlation(corr_forest, histogram = TRUE, method = "spearman")

# Let's try with the log-transformed
corr_forest_transformed <- sites %>% 
  filter(Location == "Urban Forest") %>% 
  dplyr::select(tree100m, imperv100m, log_TotalSub, 
         log_tdr, log_tdh, 
         spec_rad, light_rad, log_patch, road_length_m, trail_length_m)

findCorrelation(cor(corr_forest_transformed, method = "spearman"), cutoff = .6, verbose = TRUE, names = TRUE) 
## Compare row 10  and column  2 with corr  0.696 
##   Means:  0.443 vs 0.315 so flagging column 10 
## Compare row 2  and column  9 with corr  0.899 
##   Means:  0.361 vs 0.286 so flagging column 2 
## Compare row 9  and column  1 with corr  0.679 
##   Means:  0.279 vs 0.272 so flagging column 9 
## Compare row 6  and column  7 with corr  0.685 
##   Means:  0.371 vs 0.258 so flagging column 6 
## Compare row 7  and column  3 with corr  0.911 
##   Means:  0.296 vs 0.219 so flagging column 7 
## All correlations <= 0.6
## [1] "trail_length_m" "imperv100m"     "road_length_m"  "spec_rad"      
## [5] "light_rad"
# The same variables are dropped and kept

corr_forest_transformed <- sites %>% 
  filter(Location == "Urban Forest") %>% 
  dplyr::select(tree100m, log_TotalSub, log_tdr, patch_area_km, log_tdh)
chart.Correlation(corr_forest_transformed, histogram = TRUE, method = "spearman")

corr_center <- sites %>% 
  filter(Location == "Urban Center") %>% 
  dplyr::select(tree100m, imperv100m, TotalSub, Traffic_Dist_Road, Traffic_Dist_Highway, spec_rad, light_rad, patch_area_km, road_length_m)
# Trail length is removed because we could not measure trail length on campus

findCorrelation(cor(corr_center, method = "spearman"), cutoff = .6, verbose = TRUE, names = TRUE) 
## Compare row 1  and column  2 with corr  0.754 
##   Means:  0.477 vs 0.322 so flagging column 1 
## Compare row 2  and column  4 with corr  0.636 
##   Means:  0.388 vs 0.282 so flagging column 2 
## Compare row 7  and column  5 with corr  0.951 
##   Means:  0.352 vs 0.252 so flagging column 7 
## All correlations <= 0.6
## [1] "tree100m"   "imperv100m" "light_rad"
# This suggests that we remove tree100m, imperv100m, and light_rad
# We will keep TotalSub, Traffic_Dist_Road, Traffic_Dist_Highway, spec_rad, patch_area_km, and road_length_m

chart.Correlation(corr_center, histogram = TRUE, method = "spearman")

# Let's try with the transformed
corr_center_transformed <- sites %>% 
  filter(Location == "Urban Center") %>% 
  dplyr::select(tree100m, imperv100m, log_TotalSub, log_tdr, log_tdh, spec_rad, light_rad, log_patch, road_length_m)
findCorrelation(cor(corr_center_transformed, method = "spearman"), cutoff = .6, verbose = TRUE, names = TRUE) 
## Compare row 1  and column  2 with corr  0.754 
##   Means:  0.477 vs 0.322 so flagging column 1 
## Compare row 2  and column  4 with corr  0.636 
##   Means:  0.388 vs 0.282 so flagging column 2 
## Compare row 7  and column  5 with corr  0.951 
##   Means:  0.352 vs 0.252 so flagging column 7 
## All correlations <= 0.6
## [1] "tree100m"   "imperv100m" "light_rad"
chart.Correlation(corr_center_transformed, histogram = TRUE, method = "pearson")

# The same variables are kept or removed

For overall analysis of predictors, we will include:

For the urban forest subset, we will use:

Finally, for the urban center subset, we will use:

Graphs, Stats, and Assumptions Comparing Predictors by Location

Here, we just want to see if each predictor differs between the two Locations: urban forest and urban center. We build models, check assumptions of the model, and make graphs to represent differences.

## 
## Call:
## glm(formula = tree100m ~ Location, family = "poisson", data = sites)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4483  -1.7389  -0.8778   0.9780   2.8970  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           3.85340    0.04605   83.68   <2e-16 ***
## LocationUrban Center -3.39388    0.23399  -14.50   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 664.364  on 21  degrees of freedom
## Residual deviance:  58.556  on 20  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 6
## # A tibble: 2 × 3
##   Location      mean    se
##   <fct>        <dbl> <dbl>
## 1 Urban Forest 47.2  3.26 
## 2 Urban Center  1.58 0.692
##  Location      rate    SE  df asymp.LCL asymp.UCL
##  Urban Forest 47.15 2.171 Inf     43.08     51.61
##  Urban Center  1.58 0.363 Inf      1.01      2.48
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale

Tree cover is significantly higher in the urban forest than the urban center (z = -14.50, df = 1, 20, p < 0.001). The means and SE’s are 47.15 ± 2.17 % for the urban forest and 1.58 ± 0.36 % for the urban center.

## 
## Call:
## glm(formula = imperv100m ~ Location, family = "poisson", data = sites)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.92734  -1.59531   0.02675   0.71981   2.67358  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            0.2410     0.2803    0.86     0.39    
## LocationUrban Center   4.0965     0.2823   14.51   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1039.415  on 21  degrees of freedom
## Residual deviance:   40.746  on 20  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 6
## # A tibble: 2 × 3
##   Location      mean    se
##   <fct>        <dbl> <dbl>
## 1 Urban Forest  1.27 0.644
## 2 Urban Center 76.5  2.70
##  Location      rate    SE  df asymp.LCL asymp.UCL
##  Urban Forest  1.27 0.357 Inf     0.735       2.2
##  Urban Center 76.52 2.525 Inf    71.725      81.6
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale

Impervious cover is significantly higher in the urban center than the urban forest (z = 14.51, df = 1, 20, p < 0.001). The means and SE’s are 1.27 ± 0.36 % for the urban forest and 76.52 ± 2.53 % for the urban center.

## 
## Call:
## glm(formula = log_tdr ~ Location, family = "poisson", data = sites)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.07303  -0.43708  -0.08942   0.38719   1.04170  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            1.3786     0.1587   8.685   <2e-16 ***
## LocationUrban Center   0.2841     0.2025   1.403    0.161    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 8.5848  on 21  degrees of freedom
## Residual deviance: 6.5798  on 20  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 4
## # A tibble: 2 × 3
##   Location      mean    se
##   <fct>        <dbl> <dbl>
## 1 Urban Forest  3.97 0.349
## 2 Urban Center  5.27 0.393
##  Location     rate    SE  df asymp.LCL asymp.UCL
##  Urban Forest 3.97 0.630 Inf      2.91      5.42
##  Urban Center 5.27 0.663 Inf      4.12      6.75
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale

There is no difference in road disturbance between the urban forest and urban center (z = 1.403, df = 1, 20, p = 0.161). The means and SE’s are 52.9845308 ± 1.8776106 vehicles/day/m for the urban forest and 194.4159624 ± 1.9406054 vehicles/day/m for the urban center.

## 
## Call:
## glm(formula = log_tdh ~ Location, family = "poisson", data = sites)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.41615  -0.19952  -0.00187   0.12375   0.69255  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            1.2585     0.1685   7.467 8.22e-14 ***
## LocationUrban Center   0.1178     0.2224   0.530    0.596    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1.6528  on 21  degrees of freedom
## Residual deviance: 1.3707  on 20  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 4
## # A tibble: 2 × 3
##   Location      mean    se
##   <fct>        <dbl> <dbl>
## 1 Urban Forest  3.52 0.116
## 2 Urban Center  3.96 0.179
##  Location     rate    SE  df asymp.LCL asymp.UCL
##  Urban Forest 3.52 0.593 Inf      2.53      4.90
##  Urban Center 3.96 0.574 Inf      2.98      5.26
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale

There is no difference in highway disturbance between the urban forest and urban center (z = 0.530, df = 1, 20, p = 0.596). The means and SE’s are 33.7844285 ± 1.8094085 vehicles/day/m for the urban forest and 52.4573259 ± 1.7753543 vehicles/day/m for the urban center.

## 
## Call:
## glm(formula = log_TotalSub ~ Location, family = "poisson", data = sites)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.47861  -0.36634   0.02223   0.26908   0.61117  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            0.9173     0.1999   4.589 4.45e-06 ***
## LocationUrban Center  -0.8283     0.3409  -2.430   0.0151 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 10.6306  on 21  degrees of freedom
## Residual deviance:  4.3622  on 20  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 4
## # A tibble: 2 × 3
##   Location      mean    se
##   <fct>        <dbl> <dbl>
## 1 Urban Forest  2.50 0.139
## 2 Urban Center  1.09 0.148
##  Location     rate    SE  df asymp.LCL asymp.UCL
##  Urban Forest 2.50 0.500 Inf     1.691      3.70
##  Urban Center 1.09 0.302 Inf     0.636      1.88
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale

Plant species richness is significantly higher in the urban forest than the urban center (z = -2.430, df = 1, 20, p = 0.015). The means and SE’s are 12.182494 ± 1.6487213 species for the urban forest and 2.9742741 ± 1.3525612 species for the urban center.

## 
## Call:
## glm(formula = spec_rad ~ Location, family = "poisson", data = sites)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.73540  -0.10657  -0.01758   0.09970   0.75755  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           4.92321    0.02697 182.520  < 2e-16 ***
## LocationUrban Center  0.10975    0.03565   3.079  0.00208 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 11.0025  on 21  degrees of freedom
## Residual deviance:  1.4789  on 20  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 3
## # A tibble: 2 × 3
##   Location      mean    se
##   <fct>        <dbl> <dbl>
## 1 Urban Forest  137. 0.434
## 2 Urban Center  153. 1.26
##  Location     rate   SE  df asymp.LCL asymp.UCL
##  Urban Forest  137 3.71 Inf       130       145
##  Urban Center  153 3.58 Inf       147       161
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale

Spectral radiance is significantly higher in the urban center than the urban forest (z = 3.079, df = 1, 20, p = 0.002). The means and SE’s are 137.44 ± 3.71 Watts/(m² * sr * µm) for the urban forest and 153.39 ± 3.58 Watts/(m² * sr * µm) for the urban center.

## 
## Call:
## glm(formula = light_rad ~ Location, family = "poisson", data = sites)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.1501  -1.0590  -0.5849   1.8735   4.9674  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            1.7725     0.1303   13.60   <2e-16 ***
## LocationUrban Center   2.9888     0.1331   22.46   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1421.64  on 21  degrees of freedom
## Residual deviance:  121.73  on 20  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 4
## # A tibble: 2 × 3
##   Location       mean    se
##   <fct>         <dbl> <dbl>
## 1 Urban Forest   5.89 0.956
## 2 Urban Center 117.   9.80
##  Location       rate    SE  df asymp.LCL asymp.UCL
##  Urban Forest   5.89 0.767 Inf      4.56       7.6
##  Urban Center 116.90 3.121 Inf    110.94     123.2
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale

Light radiance is significantly higher in the urban center than the urban forest (z = 22.46, df = 1, 20, p < 0.001). The means and SE’s are 5.89 ± 0.77 mcd/m² for the urban forest and 116.90 ± 3.12 mcd/m² for the urban center.

## 
## Call:
## glm(formula = log_patch ~ Location, family = "poisson", data = sites)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.74174  -0.19285   0.01795   0.11564   0.56118  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           2.35374    0.09747  24.147  < 2e-16 ***
## LocationUrban Center -0.54815    0.15231  -3.599  0.00032 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 15.4711  on 21  degrees of freedom
## Residual deviance:  2.2471  on 20  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 4
## # A tibble: 2 × 3
##   Location      mean    se
##   <fct>        <dbl> <dbl>
## 1 Urban Forest 10.5  0.224
## 2 Urban Center  6.08 0.286
##  Location      rate    SE  df asymp.LCL asymp.UCL
##  Urban Forest 10.52 1.026 Inf      8.69     12.74
##  Urban Center  6.08 0.712 Inf      4.84      7.65
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale

Patch area is significantly higher in the urban forest than the urban center (z = -3.599, df = 1, 20, p < 0.001). The means and SE’s are 3.7049124^{4} ± 2.789884 mm² for the urban forest and 437.0291947 ± 2.0380633 mm² for the urban center.

## 
## Call:
## glm(formula = road_length_m ~ Location, family = "poisson", data = sites)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -5.585  -5.585  -1.065   2.315   8.101  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           2.74702    0.08007   34.31   <2e-16 ***
## LocationUrban Center  1.78585    0.08548   20.89   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1098.64  on 21  degrees of freedom
## Residual deviance:  446.12  on 20  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 6
## # A tibble: 2 × 3
##   Location      mean    se
##   <fct>        <dbl> <dbl>
## 1 Urban Forest  15.6  7.97
## 2 Urban Center  93.0  6.93
##  Location     rate   SE  df asymp.LCL asymp.UCL
##  Urban Forest 15.6 1.25 Inf      13.3      18.2
##  Urban Center 93.0 2.78 Inf      87.7      98.6
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale

Total road length is significantly higher in the urban center than the urban forest (z = 20.89, df = 1, 20, p < 0.001). The means and SE’s are 15.67 ± 1.25 m for the urban forest and 93.03 ± 2.78 m for the urban center.

Response Variable: Search Distance

Search Distance by Location

For this variable, I’ve gone back and forth on whether I should log-transform the search distance (because it can be ran Gaussian), or whether to use a log-link function (Poisson Distribution). The poisson glm more closely matches the mean, but the log-transformation is better at capturing the variance. Each step of the way, I have performed each test to try and gauge a best route forward. The poisson glm suggests a significant difference while the log-transformed shows no difference.

ggarrange(ggplot(sites, aes(x = WalkDist)) +
            geom_density(),
          ggplot(sites, aes(x = log(WalkDist))) +
            geom_density(), ncol = 2)

shapiro.test(sites$WalkDist) # not normal
## 
##  Shapiro-Wilk normality test
## 
## data:  sites$WalkDist
## W = 0.86946, p-value = 0.007649
shapiro.test(log(sites$WalkDist)) # normal
## 
##  Shapiro-Wilk normality test
## 
## data:  log(sites$WalkDist)
## W = 0.94443, p-value = 0.244
search <- glm(WalkDist ~ Location, data = sites, family = poisson)
summary(search)
## 
## Call:
## glm(formula = WalkDist ~ Location, family = poisson, data = sites)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -14.092   -7.441   -2.375    3.649   12.920  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           4.69410    0.03025   155.2   <2e-16 ***
## LocationUrban Center -0.46605    0.04615   -10.1   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1489.7  on 21  degrees of freedom
## Residual deviance: 1386.4  on 20  degrees of freedom
## AIC: 1517.2
## 
## Number of Fisher Scoring iterations: 5
search2 <- glm(log(WalkDist) ~ Location, data = sites)
summary(search2)
## 
## Call:
## glm(formula = log(WalkDist) ~ Location, data = sites)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.2849  -0.6105   0.2763   0.9588   1.6204  
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            3.9780     0.4062   9.792  4.5e-09 ***
## LocationUrban Center  -0.1102     0.5501  -0.200    0.843    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 1.650385)
## 
##     Null deviance: 33.074  on 21  degrees of freedom
## Residual deviance: 33.008  on 20  degrees of freedom
## AIC: 77.359
## 
## Number of Fisher Scoring iterations: 2
sites %>% 
  group_by(Location) %>% 
  summarize(mean = mean(WalkDist),
            se = plotrix::std.error(WalkDist))
## # A tibble: 2 × 3
##   Location      mean    se
##   <fct>        <dbl> <dbl>
## 1 Urban Forest 109.   31.7
## 2 Urban Center  68.6  16.2

Using a Poisson distribution, we searched significantly further at sites in the urban forest than the urban center to reach the focal web (z = -10.1, df = 1, 20, p < 0.001).

Search Distance by Predictors

# This is one way to look at the overall data, but I think I end up not using this information
models = list(WalkDist ~ 1,
              WalkDist ~ log_TotalSub, 
              WalkDist ~ log_tdr, 
              WalkDist ~ log_tdh,
              WalkDist ~ log_TotalSub + log_tdr,
              WalkDist ~ log_TotalSub + log_tdh,
              WalkDist ~ log_tdr + log_tdh,
              WalkDist ~ log_TotalSub + log_tdr + log_tdh)
fits = lapply(models, glm, data = sites, family = "poisson")
modnames = sapply(models, function(ff)deparse(ff[[3]]))
pander(aictab(fits, modname = modnames), caption="Table 1. AICc model selection table for Search Distance to the Focal Web Overall.", split.tables = Inf)
Table 1. AICc model selection table for Search Distance to the Focal Web Overall.
  Modnames K AICc Delta_AICc ModelLik AICcWt LL Cum.Wt
8 log_TotalSub + log_tdr + log_tdh 4 1504 0 1 0.9999 -747.1 0.9999
5 log_TotalSub + log_tdr 3 1523 18.94 7.709e-05 7.709e-05 -758.1 1
6 log_TotalSub + log_tdh 3 1558 53.06 3.013e-12 3.013e-12 -775.1 1
2 log_TotalSub 2 1565 60.03 9.222e-14 9.221e-14 -779.9 1
3 log_tdr 2 1597 92.5 8.208e-21 8.207e-21 -796.2 1
7 log_tdr + log_tdh 3 1600 95.14 2.187e-21 2.187e-21 -796.2 1
1 1 1 1619 114.2 1.616e-25 1.616e-25 -808.2 1
4 log_tdh 2 1621 116.6 4.912e-26 4.911e-26 -808.2 1
# The top model keeps all 3 variables and no other model is within 2 Model Likelihood Units
summary(glm(WalkDist ~ log_TotalSub + log_tdr + log_tdh, data = sites, family = poisson))
## 
## Call:
## glm(formula = WalkDist ~ log_TotalSub + log_tdr + log_tdh, family = poisson, 
##     data = sites)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -14.463   -7.694   -1.934    5.088   15.497  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   2.45699    0.24451  10.049  < 2e-16 ***
## log_TotalSub  0.31671    0.03252   9.740  < 2e-16 ***
## log_tdr       0.13279    0.01775   7.480 7.42e-14 ***
## log_tdh       0.21434    0.04469   4.796 1.62e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1489.7  on 21  degrees of freedom
## Residual deviance: 1367.3  on 18  degrees of freedom
## AIC: 1502.1
## 
## Number of Fisher Scoring iterations: 5
global <- glm(WalkDist ~ log_TotalSub + log_tdr + log_tdh, data = sites, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_sd_o <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_sd_o)
out
## 
## Call:
## glm(formula = WalkDist ~ log_tdh + log_tdr + log_TotalSub + 1, 
##     family = poisson, data = sites)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -14.463   -7.694   -1.934    5.088   15.497  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   2.45699    0.24451  10.049  < 2e-16 ***
## log_tdh       0.21434    0.04469   4.796 1.62e-06 ***
## log_tdr       0.13279    0.01775   7.480 7.42e-14 ***
## log_TotalSub  0.31671    0.03252   9.740  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1489.7  on 21  degrees of freedom
## Residual deviance: 1367.3  on 18  degrees of freedom
## AIC: 1502.1
## 
## Number of Fisher Scoring iterations: 5
info <- out[[12]]
#TotalSub (log_TotalSub)
exp(info[4, 1])
## [1] 1.372605
exp(info[4, 2])
## [1] 1.033052
#Traffic_Dist_Road (log_tdr)
exp(info[3, 1])
## [1] 1.142011
exp(info[3, 2])
## [1] 1.017911
#Traffic_Dist_Highway (log_tdh)
exp(info[2, 1])
## [1] 1.239044
exp(info[2, 2])
## [1] 1.045706
with(summary(top_model_sd_o), 1 - deviance/null.deviance) # This gives the R squared value
## [1] 0.0821148
#summary(model.avg(model_dredge, subset = delta <= 2))
#summary(model.avg(model_dredge))
model.sel(model_dredge) #estimates same signs
## Global model call: glm(formula = WalkDist ~ log_TotalSub + log_tdr + log_tdh, family = poisson, 
##     data = sites)
## ---
## Model selection table 
##   (Int)   log_tdh log_tdr log_TtS df   logLik   AICc  delta weight
## 8 2.457  0.214300 0.13280  0.3167  4 -747.072 1504.5   0.00      1
## 7 3.489           0.11210  0.2462  3 -758.052 1523.4  18.94      0
## 6 3.488  0.142400          0.2455  3 -775.110 1557.6  53.06      0
## 5 4.094                    0.2063  2 -779.947 1564.5  60.03      0
## 3 4.079           0.08155          2 -796.182 1597.0  92.50      0
## 4 4.040  0.009906 0.08199          3 -796.153 1599.6  95.14      0
## 1 4.467                            1 -808.233 1618.7 114.17      0
## 2 4.502 -0.009444                  2 -808.208 1621.0 116.55      0
## Models ranked by AICc(x)
car::vif(get.models(model_dredge, 1)[[1]]) #looking for less than 2
##      log_tdh      log_tdr log_TotalSub 
##     1.301594     1.124043     1.357283
# urban forest
global <- glm(WalkDist ~ tree100m + log_TotalSub + log_tdh + log_patch + log_tdr, data = sites_forest, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_sd_f <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_sd_f)
out
## 
## Call:
## glm(formula = WalkDist ~ log_tdh + log_tdr + log_TotalSub + tree100m + 
##     1, family = poisson, data = sites_forest)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7        8  
## -6.2159   5.9957  -8.6493   8.7999   4.4948  -2.1821  -0.1235  -9.7409  
##       9       10  
## -6.4573   5.3475  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.962190   0.608574  -4.867 1.13e-06 ***
## log_tdh       1.354582   0.109952  12.320  < 2e-16 ***
## log_tdr       0.718623   0.038386  18.721  < 2e-16 ***
## log_TotalSub -1.010040   0.095179 -10.612  < 2e-16 ***
## tree100m      0.049914   0.004241  11.770  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 892.72  on 9  degrees of freedom
## Residual deviance: 416.99  on 5  degrees of freedom
## AIC: 485.27
## 
## Number of Fisher Scoring iterations: 5
info <- out[[12]]
#tree100m
info[5, 1]
## [1] 0.04991404
info[5, 2]
## [1] 0.004240799
#TotalSub
-exp(info[4, 1])
## [1] -0.3642043
exp(info[4, 2])
## [1] 1.099856
#Traffic_Dist_Road
exp(info[3, 1])
## [1] 2.051607
exp(info[3, 2])
## [1] 1.039132
#Traffic_Dist_Highway
exp(info[2, 1])
## [1] 3.875139
exp(info[2, 2])
## [1] 1.116225
with(summary(top_model_sd_f), 1 - deviance/null.deviance) # This gives the R squared value
## [1] 0.5328974
model.sel(model_dredge) #estimates same signs
## Global model call: glm(formula = WalkDist ~ tree100m + log_TotalSub + log_tdh + 
##     log_patch + log_tdr, family = poisson, data = sites_forest)
## ---
## Model selection table 
##      (Int)   log_ptc log_tdh log_tdr log_TtS      t10 df   logLik  AICc  delta
## 31 -2.9620            1.3550  0.7186 -1.0100 0.049910  5 -237.637 500.3   0.00
## 32 -4.2870  0.111800  1.3660  0.7309 -0.9987 0.050600  6 -235.218 510.4  10.16
## 23 -4.9080            1.4220  0.5793         0.045490  4 -298.101 612.2 111.93
## 24 -6.3640  0.121500  1.4390  0.5930         0.046770  5 -295.358 615.7 115.44
## 15  1.3710            0.9436  0.4754 -0.8079           4 -311.229 638.5 138.18
## 16 -0.1680  0.130700  0.9870  0.4797 -0.8106           5 -308.126 641.3 140.98
## 29  3.7220                    0.6233 -1.1830 0.027800  4 -332.487 681.0 180.70
## 30  3.4360  0.025250          0.6257 -1.1830 0.028000  5 -332.338 689.7 189.40
## 7  -0.7565            1.1140  0.3608                   3 -355.948 721.9 221.62
## 8  -2.0680  0.111200  1.1500  0.3635                   4 -353.705 723.4 223.14
## 13  5.3650                    0.4481 -1.0220           3 -368.601 747.2 246.93
## 14  5.1990  0.015970          0.4485 -1.0240           4 -368.545 753.1 252.82
## 21  2.0920                    0.3613         0.023610  3 -413.559 837.1 336.84
## 22  1.9760  0.010040          0.3621         0.023740  4 -413.537 843.1 342.80
## 27  1.6500            0.9880         -0.3968 0.010370  4 -414.777 845.6 345.28
## 11  2.4340            0.8922         -0.3759           3 -418.557 847.1 346.84
## 12  2.0780  0.029380  0.9042         -0.3742           4 -418.359 852.7 352.44
## 28  1.7170 -0.006520  0.9872         -0.3978 0.010520  5 -414.767 854.5 354.26
## 19  0.7985            0.9742                 0.008801  3 -428.397 866.8 366.52
## 3   1.4210            0.9170                           2 -431.230 868.2 367.90
## 4   1.0170  0.034130  0.9295                           3 -430.969 871.9 371.66
## 20  0.7048  0.009042  0.9759                 0.008645  4 -428.379 872.8 372.48
## 5   3.7340                    0.2338                   2 -437.478 880.7 380.40
## 6   3.9550 -0.021130          0.2341                   3 -437.380 884.8 384.49
## 9   5.7640                           -0.4343           2 -457.538 920.8 420.52
## 10  6.1830 -0.039100                 -0.4371           3 -457.142 924.3 424.01
## 25  5.6580                           -0.4322 0.002144  3 -457.300 924.6 424.33
## 26  6.1230 -0.046060                 -0.4359 0.002736  4 -456.763 929.5 429.25
## 1   4.6940                                             1 -475.501 953.5 453.23
## 17  4.5250                                   0.003571  2 -474.832 955.4 455.10
## 2   5.0660 -0.035390                                   2 -475.193 956.1 455.83
## 18  4.9430 -0.041030                         0.003855  3 -474.415 958.8 458.56
##    weight
## 31  0.994
## 32  0.006
## 23  0.000
## 24  0.000
## 15  0.000
## 16  0.000
## 29  0.000
## 30  0.000
## 7   0.000
## 8   0.000
## 13  0.000
## 14  0.000
## 21  0.000
## 22  0.000
## 27  0.000
## 11  0.000
## 12  0.000
## 28  0.000
## 19  0.000
## 3   0.000
## 4   0.000
## 20  0.000
## 5   0.000
## 6   0.000
## 9   0.000
## 10  0.000
## 25  0.000
## 26  0.000
## 1   0.000
## 17  0.000
## 2   0.000
## 18  0.000
## Models ranked by AICc(x)
#summary(model.avg(model_dredge, subset = delta <= 2))
car::vif(get.models(model_dredge, 1)[[1]]) #looking for less than 2
##      log_tdh      log_tdr log_TotalSub     tree100m 
##     1.444130     1.612738     1.259335     1.840585
# urban center
global <- glm(WalkDist ~ log_TotalSub + log_tdr + log_tdh + spec_rad + log_patch + road_length_m, data = sites_center, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_sd_c <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_sd_c)
out
## 
## Call:
## glm(formula = WalkDist ~ log_patch + log_tdh + log_tdr + log_TotalSub + 
##     spec_rad + 1, family = poisson, data = sites_center)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3145  -1.2020  -0.2198   0.6284   3.2692  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -26.90656    2.78649  -9.656  < 2e-16 ***
## log_patch      1.05254    0.05903  17.832  < 2e-16 ***
## log_tdh        0.58001    0.09128   6.354 2.10e-10 ***
## log_tdr        0.45331    0.04386  10.335  < 2e-16 ***
## log_TotalSub   0.54581    0.08135   6.709 1.96e-11 ***
## spec_rad       0.12466    0.01458   8.549  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 493.631  on 11  degrees of freedom
## Residual deviance:  29.624  on  6  degrees of freedom
## AIC: 110.15
## 
## Number of Fisher Scoring iterations: 4
info <- out[[12]]
#TotalSub
exp(info[5, 1])
## [1] 1.726
exp(info[5, 2])
## [1] 1.084755
#patch_area_km
exp(info[2, 1])
## [1] 2.864928
exp(info[2, 2])
## [1] 1.060802
#spec_rad
info[6, 1]
## [1] 0.1246571
info[6, 2]
## [1] 0.01458074
#Traffic_Dist_Road
exp(info[4, 1])
## [1] 1.573514
exp(info[4, 2])
## [1] 1.044837
#Traffic_Dist_Highway
exp(info[3, 1])
## [1] 1.786049
exp(info[3, 2])
## [1] 1.095576
with(summary(top_model_sd_c), 1 - deviance/null.deviance) # This gives the R squared value
## [1] 0.9399876
model.sel(model_dredge) #estimates same signs
## Global model call: glm(formula = WalkDist ~ log_TotalSub + log_tdr + log_tdh + spec_rad + 
##     log_patch + road_length_m, family = poisson, data = sites_center)
## ---
## Model selection table 
##        (Int) log_ptc  log_tdh log_tdr  log_TtS  rod_lng_m   spc_rad df   logLik
## 48 -26.91000  1.0530  0.58000  0.4533 0.545800             0.124700  6  -49.077
## 64 -27.07000  1.0540  0.58820  0.4582 0.544400  4.438e-04  0.125100  7  -49.045
## 46 -18.79000  0.9700           0.3030 0.467300             0.096030  5  -69.398
## 62 -18.05000  0.9688           0.2870 0.477700 -2.947e-03  0.093480  6  -67.840
## 40 -23.31000  0.9764  0.50860  0.5486                      0.106700  5  -72.569
## 56 -23.55000  0.9785  0.53010  0.5595           9.858e-04  0.106600  6  -72.413
## 38 -16.35000  0.9141           0.4031                      0.082200  4  -87.897
## 54 -16.16000  0.9160           0.3926          -2.429e-03  0.082660  5  -86.790
## 16  -4.79700  0.8060  0.35770  0.3880 0.347000                       5  -92.534
## 32  -4.18200  0.7949  0.30640  0.3545 0.376000 -2.056e-03            6  -91.779
## 30  -2.16500  0.7837           0.2588 0.371000 -4.289e-03            5  -97.463
## 14  -2.78000  0.7983           0.2954 0.306900                       4 -101.490
## 8   -4.47000  0.7813  0.32270  0.4521                                4 -102.558
## 24  -4.43700  0.7806  0.31990  0.4506          -9.962e-05            5 -102.556
## 6   -2.69600  0.7794           0.3634                                3 -109.821
## 58 -11.68000  0.8625                  0.840500 -6.534e-03  0.066040  5 -104.636
## 22  -2.30800  0.7684           0.3491          -2.607e-03            4 -108.248
## 60 -11.77000  0.8623  0.01075         0.845300 -6.509e-03  0.066310  6 -104.625
## 42 -13.35000  0.8539                  0.825900             0.073500  4 -114.026
## 44 -13.62000  0.8523  0.03948         0.841100             0.074180  5 -113.881
## 26  -0.33890  0.7077                  0.716900 -7.737e-03            4 -124.349
## 28  -0.04446  0.7153 -0.07423         0.688600 -7.937e-03            5 -123.800
## 10  -0.84410  0.6920                  0.625600                       3 -140.065
## 12  -0.73870  0.6951 -0.02912         0.616400                       4 -139.983
## 52  -2.37500  0.6594 -0.24640                  -6.566e-03  0.026110  5 -163.515
## 20   1.97600  0.6132 -0.25760                  -6.713e-03            4 -166.678
## 50  -3.52000  0.6142                           -5.326e-03  0.028380  4 -169.190
## 18   1.17400  0.5606                           -5.367e-03            3 -173.031
## 36  -3.50600  0.6315 -0.17530                              0.028930  4 -173.540
## 34  -4.11800  0.6026                                       0.029610  3 -176.457
## 4    1.21900  0.5873 -0.17900                                        3 -177.822
## 2    0.72240  0.5547                                                 2 -180.946
## 21   3.78000                   0.1896          -6.382e-03            3 -246.767
## 53   6.51000                   0.1793          -6.883e-03 -0.017160  4 -245.220
## 29   3.80300                   0.1709 0.122100 -7.024e-03            4 -245.751
## 23   3.34400          0.07582  0.2052          -5.828e-03            4 -246.089
## 61   6.43400                   0.1612 0.117700 -7.549e-03 -0.016480  5 -244.290
## 55   5.95200          0.06141  0.1931          -6.387e-03 -0.015890  5 -244.790
## 31   3.42200          0.06619  0.1858 0.113100 -6.496e-03            5 -245.228
## 7    2.43700          0.14180  0.2262                                3 -252.264
## 5    3.15300                   0.1976                                2 -254.778
## 63   5.96200          0.05205  0.1739 0.110600 -7.086e-03 -0.015420  6 -243.978
## 39   3.39400          0.14100  0.2231                     -0.006122  4 -252.053
## 37   4.17700                   0.1940                     -0.006552  3 -254.531
## 15   2.43500          0.14110  0.2235 0.017190                       4 -252.240
## 13   3.14700                   0.1940 0.023060                       3 -254.735
## 45   4.12500                   0.1922 0.013050            -0.006248  4 -254.518
## 47   3.36500          0.14070  0.2218 0.008304            -0.005939  5 -252.048
## 57   8.67200                          0.317300 -9.203e-03 -0.025840  4 -257.913
## 59   9.27400         -0.09140         0.304500 -9.818e-03 -0.026950  5 -256.705
## 25   4.59900                          0.343500 -8.289e-03            3 -262.339
## 27   4.96100         -0.07796         0.334600 -8.768e-03            4 -261.447
## 49   9.56200                                   -7.437e-03 -0.030380  3 -265.893
## 51  10.32000         -0.11640                  -8.312e-03 -0.031830  4 -264.034
## 17   4.82400                                   -6.523e-03            2 -271.928
## 19   5.27900         -0.09958                  -7.195e-03            3 -270.558
## 9    3.96600                          0.233400                       2 -275.865
## 41   6.01800                          0.202000            -0.013160  3 -274.663
## 11   4.02500         -0.01410         0.231200                       3 -275.836
## 33   7.22500                                              -0.019560  2 -278.356
## 43   6.06900         -0.01300         0.200000            -0.013140  4 -274.639
## 1    4.22800                                                         1 -281.081
## 35   7.33500         -0.03213                             -0.019450  3 -278.209
## 3    4.36800         -0.03536                                        2 -280.903
##     AICc  delta weight
## 48 127.0   0.00  0.999
## 64 140.1  13.14  0.001
## 46 158.8  31.84  0.000
## 62 164.5  37.53  0.000
## 40 165.1  38.18  0.000
## 56 173.6  46.67  0.000
## 38 189.5  62.55  0.000
## 54 193.6  66.63  0.000
## 16 205.1  78.11  0.000
## 32 212.4  85.40  0.000
## 30 214.9  87.97  0.000
## 14 216.7  89.74  0.000
## 8  218.8  91.88  0.000
## 24 225.1  98.16  0.000
## 6  228.6 101.69  0.000
## 58 229.3 102.32  0.000
## 22 230.2 103.26  0.000
## 60 238.0 111.10  0.000
## 42 241.8 114.81  0.000
## 44 247.8 120.81  0.000
## 26 262.4 135.46  0.000
## 28 267.6 140.65  0.000
## 10 289.1 162.18  0.000
## 12 293.7 166.73  0.000
## 52 347.0 220.08  0.000
## 20 347.1 220.12  0.000
## 50 352.1 225.14  0.000
## 18 355.1 228.11  0.000
## 36 360.8 233.84  0.000
## 34 361.9 234.96  0.000
## 4  364.6 237.69  0.000
## 2  367.2 240.27  0.000
## 21 502.5 375.58  0.000
## 53 504.2 377.20  0.000
## 29 505.2 378.26  0.000
## 23 505.9 378.94  0.000
## 61 508.6 381.63  0.000
## 55 509.6 382.63  0.000
## 31 510.5 383.50  0.000
## 7  513.5 386.57  0.000
## 5  514.9 387.94  0.000
## 63 516.8 389.80  0.000
## 39 517.8 390.87  0.000
## 37 518.1 391.11  0.000
## 15 518.2 391.24  0.000
## 13 518.5 391.52  0.000
## 45 522.8 395.80  0.000
## 47 524.1 397.14  0.000
## 57 529.5 402.59  0.000
## 59 533.4 406.46  0.000
## 25 533.7 406.72  0.000
## 27 536.6 409.65  0.000
## 49 540.8 413.83  0.000
## 51 541.8 414.83  0.000
## 17 549.2 422.24  0.000
## 19 550.1 423.16  0.000
## 9  557.1 430.11  0.000
## 41 558.3 431.37  0.000
## 11 560.7 433.72  0.000
## 33 562.0 435.09  0.000
## 43 563.0 436.04  0.000
## 1  564.6 437.61  0.000
## 35 565.4 438.46  0.000
## 3  567.1 440.19  0.000
## Models ranked by AICc(x)
#summary(model.avg(model_dredge, subset = delta <= 2))
car::vif(get.models(model_dredge, 1)[[1]]) # road and highway disturbance between 2 and 3
##    log_patch      log_tdh      log_tdr log_TotalSub     spec_rad 
##     1.874528     2.076776     2.677800     1.554569     1.661330

Overall Coefficients: Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.45699 0.24451 10.049 < 2e-16 log_tdh 0.21434 0.04469 4.796 1.62e-06 positive log_tdr 0.13279 0.01775 7.480 7.42e-14 *** positive log_TotalSub 0.31671 0.03252 9.740 < 2e-16 *** positive

R-squared = 0.0821148

Urban Forest Coefficients: Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.962190 0.608574 -4.867 1.13e-06 log_tdh 1.354582 0.109952 12.320 < 2e-16 positive log_tdr 0.718623 0.038386 18.721 < 2e-16 *** positive log_TotalSub -1.010040 0.095179 -10.612 < 2e-16 *** negative tree100m 0.049914 0.004241 11.770 < 2e-16 *** positive

R-squared = 0.5328974

Urban Center Coefficients: Estimate Std. Error z value Pr(>|z|)
(Intercept) -26.90656 2.78649 -9.656 < 2e-16 log_patch 1.05254 0.05903 17.832 < 2e-16 positive log_tdh 0.58001 0.09128 6.354 2.10e-10 *** positive log_tdr 0.45331 0.04386 10.335 < 2e-16 *** positive log_TotalSub 0.54581 0.08135 6.709 1.96e-11 *** positive spec_rad 0.12466 0.01458 8.549 < 2e-16 *** positive

R-squared = 0.9399876

Search Distance by Land

search_land <- glm(WalkDist ~ Land, data = sites, family = poisson)
summary(search_land)
## 
## Call:
## glm(formula = WalkDist ~ Land, family = poisson, data = sites)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -12.800   -6.567   -1.276    3.797   15.068  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           3.31419    0.09535   34.76   <2e-16 ***
## LandUrban, Medium     1.15172    0.10491   10.98   <2e-16 ***
## LandUrban, Low        1.24494    0.11969   10.40   <2e-16 ***
## LandDeciduous Forest  1.20293    0.10320   11.66   <2e-16 ***
## LandWoody Wetlands    1.70088    0.10632   16.00   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1489.7  on 21  degrees of freedom
## Residual deviance: 1149.3  on 17  degrees of freedom
## AIC: 1286.1
## 
## Number of Fisher Scoring iterations: 5
car::Anova(search_land, test.statistic = "LR")
## Analysis of Deviance Table (Type II tests)
## 
## Response: WalkDist
##      LR Chisq Df Pr(>Chisq)    
## Land   340.32  4  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
comp <- glht(search_land, linfct = mcp(Land = "Tukey"))
summary(comp)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: glm(formula = WalkDist ~ Land, family = poisson, data = sites)
## 
## Linear Hypotheses:
##                                        Estimate Std. Error z value Pr(>|z|)    
## Urban, Medium - Urban, High == 0        1.15172    0.10491  10.978   <1e-04 ***
## Urban, Low - Urban, High == 0           1.24494    0.11969  10.401   <1e-04 ***
## Deciduous Forest - Urban, High == 0     1.20293    0.10320  11.656   <1e-04 ***
## Woody Wetlands - Urban, High == 0       1.70088    0.10632  15.998   <1e-04 ***
## Urban, Low - Urban, Medium == 0         0.09322    0.08457   1.102    0.796    
## Deciduous Forest - Urban, Medium == 0   0.05121    0.05896   0.869    0.903    
## Woody Wetlands - Urban, Medium == 0     0.54916    0.06425   8.547   <1e-04 ***
## Deciduous Forest - Urban, Low == 0     -0.04201    0.08244  -0.510    0.986    
## Woody Wetlands - Urban, Low == 0        0.45594    0.08630   5.283   <1e-04 ***
## Woody Wetlands - Deciduous Forest == 0  0.49795    0.06142   8.107   <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
sites <- sites %>% 
  mutate(Land = fct_relevel(Land, "Deciduous Forest", 
                            "Woody Wetlands", 
                            "Urban, High", 
                            "Urban, Medium", 
                            "Urban, Low"))
nd <- data.frame(Land = factor(levels(sites$Land)))
predictions <- augment(search_land, newdata = nd, se_fit = TRUE, type.predict = "response")
predictions <- predictions %>% 
  rename("rate" = ".fitted", "SE" = ".se.fit")
#predictions <- summary(emmeans(search_land, ~Land),type = "response")

ggplot(sites, aes(x = Land, y = WalkDist)) + 
  geom_jitter(color = "grey", width = 0.1, size = 0.5) +
  geom_point(aes(x = Land, y = rate, color = Land), size = 1, data = predictions) + 
  geom_errorbar(aes(x = Land, 
                  ymin = rate - SE, 
                  ymax = rate + SE, color = Land), data = predictions, inherit.aes = FALSE, width = 0.25, size = 1) +
  theme_classic() +
  scale_color_manual(values = c("Deciduous Forest" = "#1b9e77", "Woody Wetlands" = "#1b9e77", "Urban, High" = "#d95f02", "Urban, Medium" = "#d95f02", "Urban, Low" = "#d95f02")) +
  xlab("Land Cover Class") +
  ylab("Search Distance [m]") + 
  theme(text = element_text(size = 10)) + 
  theme(axis.text.x=element_text(colour="black", size=10)) + 
  theme(axis.text.y=element_text(colour="black", size=10)) +
  theme(legend.position = "none") +
  theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust=1)) +
  annotate(geom="text", x=1, y=max(sites$WalkDist), label="A", color="black", size = 3) +
  annotate(geom="text", x=2, y=max(sites$WalkDist), label="B", color="black", size = 3) +
  annotate(geom="text", x=3, y=max(sites$WalkDist), label="C", color="black", size = 3) +
  annotate(geom="text", x=4, y=max(sites$WalkDist), label="A", color="black", size = 3) +
  annotate(geom="text", x=5, y=max(sites$WalkDist), label="A", color="black", size = 3)

Search distance significantly varies by land class (Chi = 340.32, df = 4, p < 0.001).

Response Variable: Number of Webs

Number of Webs by Location

I also tested Poisson glm vs log-transformation here since log-transformed does help make the data more Gaussian. However, similar results arise in this case.

ggarrange(ggplot(sites, aes(x = NumWebs)) +
            geom_density(),
          ggplot(sites, aes(x = log(NumWebs))) +
            geom_density(), ncol = 2)

shapiro.test(sites$NumWebs) # not normal
## 
##  Shapiro-Wilk normality test
## 
## data:  sites$NumWebs
## W = 0.54806, p-value = 3.839e-07
shapiro.test(log(sites$NumWebs)) # just normal
## 
##  Shapiro-Wilk normality test
## 
## data:  log(sites$NumWebs)
## W = 0.91313, p-value = 0.05487
web <- glm(NumWebs ~ Location, data = sites, family = poisson)
summary(web)
## 
## Call:
## glm(formula = NumWebs ~ Location, family = poisson, data = sites)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.9610  -2.2514  -0.3075   0.0000  10.3449  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            1.0986     0.1826   6.017 1.77e-09 ***
## LocationUrban Center   1.6792     0.1963   8.556  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 305.94  on 21  degrees of freedom
## Residual deviance: 200.80  on 20  degrees of freedom
## AIC: 283.55
## 
## Number of Fisher Scoring iterations: 5
web2 <- glm(log(NumWebs) ~ Location, data = sites)
summary(web)
## 
## Call:
## glm(formula = NumWebs ~ Location, family = poisson, data = sites)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.9610  -2.2514  -0.3075   0.0000  10.3449  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            1.0986     0.1826   6.017 1.77e-09 ***
## LocationUrban Center   1.6792     0.1963   8.556  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 305.94  on 21  degrees of freedom
## Residual deviance: 200.80  on 20  degrees of freedom
## AIC: 283.55
## 
## Number of Fisher Scoring iterations: 5
# Both methods yield similar results
sites %>% 
  group_by(Location) %>% 
  summarize(mean = mean(NumWebs),
            se = plotrix::std.error(NumWebs))
## # A tibble: 2 × 3
##   Location      mean    se
##   <fct>        <dbl> <dbl>
## 1 Urban Forest   3   0.211
## 2 Urban Center  16.1 5.76

There were significantly more webs at sites in the urban center than the urban forest (z = 8.556, df = 1, 20, p < 0.001).

Number of Webs by Predictors

# This is one way to look at the overall data, but I think I end up not using this information
models = list(NumWebs ~ 1,
              NumWebs ~ log_TotalSub, 
              NumWebs ~ log_tdr, 
              NumWebs ~ log_tdh,
              NumWebs ~ log_TotalSub + log_tdr,
              NumWebs ~ log_TotalSub + log_tdh,
              NumWebs ~ log_tdr + log_tdh,
              NumWebs ~ log_TotalSub + log_tdr + log_tdh)
fits = lapply(models, glm, data = sites, family = "poisson")
modnames = sapply(models, function(ff)deparse(ff[[3]]))
pander(aictab(fits, modname = modnames), caption="Table 1. AICc model selection table for Search Distance to the Focal Web Overall.", split.tables = Inf)
Table 1. AICc model selection table for Search Distance to the Focal Web Overall.
  Modnames K AICc Delta_AICc ModelLik AICcWt LL Cum.Wt
8 log_TotalSub + log_tdr + log_tdh 4 224.9 0 1 0.9985 -107.3 0.9985
6 log_TotalSub + log_tdh 3 237.8 12.94 0.001549 0.001546 -115.2 1
5 log_TotalSub + log_tdr 3 274.9 50 1.385e-11 1.383e-11 -133.8 1
2 log_TotalSub 2 304.1 79.17 6.434e-18 6.424e-18 -149.7 1
7 log_tdr + log_tdh 3 339.1 114.3 1.542e-25 1.539e-25 -165.9 1
3 log_tdr 2 343.7 118.8 1.563e-26 1.561e-26 -169.5 1
4 log_tdh 2 374.6 149.7 3.126e-33 3.121e-33 -185 1
1 1 1 386.9 162 6.602e-36 6.591e-36 -192.3 1
# The top model keeps all 3 variables and no other model is within 2 Model Likelihood Units
summary(glm(NumWebs ~ log_TotalSub + log_tdr + log_tdh, data = sites, family = poisson))
## 
## Call:
## glm(formula = NumWebs ~ log_TotalSub + log_tdr + log_tdh, family = poisson, 
##     data = sites)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.7503  -1.6315  -0.7943   1.0933   5.5758  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   7.55449    0.87254   8.658  < 2e-16 ***
## log_TotalSub -1.02619    0.09768 -10.505  < 2e-16 ***
## log_tdr       0.18726    0.04729   3.960 7.49e-05 ***
## log_tdh      -1.26833    0.19726  -6.430 1.28e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 305.94  on 21  degrees of freedom
## Residual deviance: 135.78  on 18  degrees of freedom
## AIC: 222.53
## 
## Number of Fisher Scoring iterations: 5
global <- glm(NumWebs ~ log_TotalSub + log_tdr + log_tdh, data = sites, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_nw_o <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_nw_o)
out
## 
## Call:
## glm(formula = NumWebs ~ log_tdh + log_tdr + log_TotalSub + 1, 
##     family = poisson, data = sites)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.7503  -1.6315  -0.7943   1.0933   5.5758  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   7.55449    0.87254   8.658  < 2e-16 ***
## log_tdh      -1.26833    0.19726  -6.430 1.28e-10 ***
## log_tdr       0.18726    0.04729   3.960 7.49e-05 ***
## log_TotalSub -1.02619    0.09768 -10.505  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 305.94  on 21  degrees of freedom
## Residual deviance: 135.78  on 18  degrees of freedom
## AIC: 222.53
## 
## Number of Fisher Scoring iterations: 5
info <- out[[12]]
#TotalSub (log_TotalSub)
exp(info[4, 1])
## [1] 0.3583704
exp(info[4, 2])
## [1] 1.102614
#Traffic_Dist_Road (log_tdr)
exp(info[3, 1])
## [1] 1.20594
exp(info[3, 2])
## [1] 1.048422
#Traffic_Dist_Highway (log_tdh)
exp(info[2, 1])
## [1] 0.2812998
exp(info[2, 2])
## [1] 1.21806
with(summary(top_model_nw_o), 1 - deviance/null.deviance) # This gives the R squared value
## [1] 0.556197
#summary(model.avg(model_dredge, subset = delta <= 2))
#summary(model.avg(model_dredge))
model.sel(model_dredge) #estimates same signs
## Global model call: glm(formula = NumWebs ~ log_TotalSub + log_tdr + log_tdh, family = poisson, 
##     data = sites)
## ---
## Model selection table 
##    (Int) log_tdh log_tdr log_TtS df   logLik  AICc  delta weight
## 8 7.5540 -1.2680  0.1873 -1.0260  4 -107.266 224.9   0.00  0.998
## 6 9.1080 -1.4250         -1.0500  3 -115.246 237.8  12.94  0.002
## 7 1.9820          0.2824 -0.7399  3 -133.778 274.9  50.00  0.000
## 5 3.4270                 -0.7551  2 -149.711 304.1  79.17  0.000
## 4 2.1490 -0.3535  0.2955          3 -165.907 339.1 114.26  0.000
## 3 0.6916          0.3255          2 -169.547 343.7 118.84  0.000
## 2 4.2430 -0.5220                  2 -184.972 374.6 149.69  0.000
## 1 2.3160                          1 -192.348 386.9 162.01  0.000
## Models ranked by AICc(x)
car::vif(get.models(model_dredge, 1)[[1]]) #looking for less than 2
##      log_tdh      log_tdr log_TotalSub 
##     1.359106     1.059895     1.291282
# urban forest
global <- glm(NumWebs ~ tree100m + log_TotalSub + log_tdh + log_patch + log_tdr, data = sites_forest, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_nw_f <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_nw_f)
out
## 
## Call:
## glm(formula = NumWebs ~ 1, family = poisson, data = sites_forest)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.6149   0.0000   0.0000   0.0000   0.5491  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.0986     0.1826   6.017 1.77e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1.3592  on 9  degrees of freedom
## Residual deviance: 1.3592  on 9  degrees of freedom
## AIC: 33.069
## 
## Number of Fisher Scoring iterations: 4
info <- out[[12]]
with(summary(top_model_nw_f), 1 - deviance/null.deviance) # This gives the R squared value
## [1] -6.661338e-16
model.sel(model_dredge) #estimates same signs
## Global model call: glm(formula = NumWebs ~ tree100m + log_TotalSub + log_tdh + log_patch + 
##     log_tdr, family = poisson, data = sites_forest)
## ---
## Model selection table 
##      (Int)  log_ptc log_tdh  log_tdr log_TtS        t10 df  logLik AICc delta
## 1   1.0990                                               1 -15.535 33.6  0.00
## 9   0.2625                            0.3303             2 -15.252 36.2  2.65
## 3   2.2630          -0.3326                              2 -15.320 36.4  2.79
## 2   0.3896 0.067270                                      2 -15.504 36.7  3.15
## 5   1.0320                   0.01662                     2 -15.530 36.8  3.21
## 17  1.1360                                   -0.0007896  2 -15.534 36.8  3.21
## 11  1.2880          -0.2564           0.2798             3 -15.125 40.2  6.68
## 13  0.3186                  -0.04425  0.3777             3 -15.224 40.4  6.88
## 10 -0.4014 0.064140                   0.3255             3 -15.226 40.5  6.88
## 25  0.1897                            0.3355  0.0012660  3 -15.250 40.5  6.93
## 7   2.6440          -0.3940 -0.04199                     3 -15.297 40.6  7.02
## 4   1.7900 0.041070 -0.3213                              3 -15.309 40.6  7.05
## 19  2.2840          -0.3319                  -0.0005167  3 -15.320 40.6  7.07
## 6   0.2612 0.071370          0.02142                     3 -15.496 41.0  7.42
## 18  0.4267 0.069970                          -0.0013900  3 -15.501 41.0  7.43
## 21  1.0200                   0.01755          0.0001927  3 -15.530 41.1  7.49
## 15  2.0160          -0.4133 -0.11500  0.3896             4 -14.981 46.0 12.39
## 12  0.7844 0.044000 -0.2446           0.2792             4 -15.112 46.2 12.66
## 27  1.1980          -0.2589           0.2881  0.0016430  4 -15.121 46.2 12.67
## 14 -0.2805 0.057190         -0.03979  0.3692             4 -15.203 46.4 12.84
## 29  0.4260                  -0.05314  0.3807 -0.0016920  4 -15.221 46.4 12.87
## 26 -0.4464 0.063020                   0.3299  0.0009705  4 -15.225 46.4 12.88
## 23  3.0710          -0.4285 -0.06889         -0.0042350  4 -15.278 46.6 12.99
## 8   2.2560 0.030320 -0.3794 -0.03762                     4 -15.291 46.6 13.01
## 20  1.8090 0.041980 -0.3198                  -0.0007209  4 -15.308 46.6 13.05
## 22  0.2718 0.071530          0.02051         -0.0001859  4 -15.496 47.0 13.42
## 31  2.6290          -0.4631 -0.15830  0.4011 -0.0062770  5 -14.943 54.9 21.32
## 16  1.8370 0.014290 -0.4063 -0.11260  0.3873             5 -14.980 55.0 21.39
## 28  0.7145 0.042870 -0.2477           0.2871  0.0015380  5 -15.109 55.2 21.65
## 30 -0.1691 0.057990         -0.04964  0.3721 -0.0018670  5 -15.200 55.4 21.83
## 24  2.6970 0.028670 -0.4133 -0.06449         -0.0041930  5 -15.272 55.5 21.98
## 32  2.5160 0.008674 -0.4581 -0.15660  0.3994 -0.0062390  6 -14.943 69.9 36.32
##    weight
## 1   0.411
## 9   0.109
## 3   0.102
## 2   0.085
## 5   0.083
## 17  0.083
## 11  0.015
## 13  0.013
## 10  0.013
## 25  0.013
## 7   0.012
## 4   0.012
## 19  0.012
## 6   0.010
## 18  0.010
## 21  0.010
## 15  0.001
## 12  0.001
## 27  0.001
## 14  0.001
## 29  0.001
## 26  0.001
## 23  0.001
## 8   0.001
## 20  0.001
## 22  0.001
## 31  0.000
## 16  0.000
## 28  0.000
## 30  0.000
## 24  0.000
## 32  0.000
## Models ranked by AICc(x)
#summary(model.avg(model_dredge, subset = delta <= 2))
#car::vif(get.models(model_dredge, 1)[[1]]) #looking for less than 2

# urban center
global <- glm(NumWebs ~ log_TotalSub + log_tdr + log_tdh + spec_rad + log_patch + road_length_m, data = sites_center, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_nw_c <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_nw_c)
out
## 
## Call:
## glm(formula = NumWebs ~ log_patch + log_tdh + log_TotalSub + 
##     road_length_m + 1, family = poisson, data = sites_center)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -3.005  -1.905  -1.369   1.040   5.011  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   13.145965   1.236110  10.635  < 2e-16 ***
## log_patch     -0.293327   0.113716  -2.579   0.0099 ** 
## log_tdh       -1.539660   0.190172  -8.096 5.67e-16 ***
## log_TotalSub  -0.717109   0.183835  -3.901 9.59e-05 ***
## road_length_m -0.022553   0.004713  -4.785 1.71e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 199.436  on 11  degrees of freedom
## Residual deviance:  83.466  on  7  degrees of freedom
## AIC: 142.51
## 
## Number of Fisher Scoring iterations: 6
info <- out[[12]]
#TotalSub
exp(info[4, 1])
## [1] 0.4881615
exp(info[4, 2])
## [1] 1.201817
#patch_area_km
exp(info[2, 1])
## [1] 0.7457781
exp(info[2, 2])
## [1] 1.120434
#Traffic_Dist_Highway
exp(info[3, 1])
## [1] 0.2144539
exp(info[3, 2])
## [1] 1.209458
#road_length
info[5, 1]
## [1] -0.02255259
info[5, 2]
## [1] 0.004713336
with(summary(top_model_nw_c), 1 - deviance/null.deviance) # This gives the R squared value
## [1] 0.5814882
model.sel(model_dredge) #estimates same signs
## Global model call: glm(formula = NumWebs ~ log_TotalSub + log_tdr + log_tdh + spec_rad + 
##     log_patch + road_length_m, family = poisson, data = sites_center)
## ---
## Model selection table 
##      (Int)  log_ptc log_tdh   log_tdr log_TtS rod_lng_m    spc_rad df   logLik
## 28 13.1500 -0.29330  -1.540           -0.7171 -0.022550             5  -66.255
## 27 11.1500           -1.589           -0.5827 -0.019480             4  -69.677
## 31 10.8400           -1.560  0.048410 -0.6566 -0.019350             5  -69.297
## 59  9.3570           -1.598           -0.5745 -0.019350  0.0117000  5  -69.534
## 32 12.8600 -0.30680  -1.499  0.059820 -0.8179 -0.022620             6  -65.643
## 60 16.0800 -0.33480  -1.515           -0.7464 -0.023170 -0.0174800  6  -65.996
## 19  9.9200           -1.404                   -0.020320             3  -76.199
## 20 11.3900 -0.20440  -1.392                   -0.023460             4  -74.349
## 23 10.4000           -1.462 -0.045980         -0.020440             4  -75.771
## 51  7.1470           -1.422                   -0.019770  0.0181500  4  -75.804
## 63  6.6550           -1.564  0.077350 -0.6834 -0.019180  0.0263300  6  -68.758
## 24 11.7900 -0.20110  -1.444 -0.042730         -0.023390             5  -73.972
## 52 11.4600 -0.20550  -1.392                   -0.023500 -0.0004335  5  -74.349
## 55  8.4240           -1.457 -0.032800         -0.020020  0.0120200  5  -75.630
## 11  9.4370           -1.598           -0.5986                       3  -81.277
## 64 13.8000 -0.31780  -1.495  0.053790 -0.8168 -0.022770 -0.0054980  7  -65.625
## 12 10.3300 -0.16630  -1.537           -0.7017                       4  -79.845
## 43  4.4710           -1.638           -0.5498            0.0329400  4  -80.113
## 15  8.9470           -1.545  0.069690 -0.6827                       4  -80.571
## 56 14.3800 -0.23710  -1.446 -0.058370         -0.024430 -0.0142400  6  -73.828
## 47  1.4920           -1.576  0.106900 -0.6659            0.0478100  5  -78.606
## 16  9.8490 -0.17370  -1.480  0.074330 -0.7962                       5  -79.029
## 44  6.8720 -0.12500  -1.582           -0.6431            0.0215200  5  -79.437
## 35  1.1250           -1.456                              0.0465400  3  -86.735
## 3   7.9570           -1.370                                         2  -89.716
## 48  3.7960 -0.11200  -1.531  0.102100 -0.7422            0.0368400  6  -78.062
## 36  0.6251  0.03279  -1.471                              0.0488700  4  -86.673
## 39  0.8578           -1.446  0.008436                    0.0477500  4  -86.724
## 7   8.3620           -1.423 -0.038240                               3  -89.465
## 4   8.1330 -0.03739  -1.356                                         3  -89.626
## 8   8.5580 -0.04012  -1.409 -0.039270                               4  -89.362
## 40  0.2035  0.03577  -1.458  0.011780                    0.0507700  5  -86.652
## 30  5.1800 -0.27760          0.215200 -0.7836 -0.011890             5 -105.578
## 29  3.0470                   0.216400 -0.5251 -0.009834             4 -109.766
## 61 -4.7150                   0.258900 -0.5583 -0.010410  0.0495400  5 -107.225
## 14  3.7230 -0.22450          0.227800 -0.7784                       4 -111.199
## 45 -5.9680                   0.269300 -0.5592            0.0512800  4 -111.473
## 13  2.1090                   0.236300 -0.5758                       3 -114.223
## 26  5.9330 -0.25700                   -0.5323 -0.011640             4 -112.193
## 21  3.0990                   0.135900         -0.011710             3 -114.554
## 53 -3.2890                   0.171400         -0.011440  0.0401500  4 -112.589
## 62  1.1030 -0.22690          0.237100 -0.7515 -0.011810  0.0235000  6 -105.169
## 22  3.9920 -0.12100          0.126900         -0.013000             4 -113.487
## 17  3.8080                                    -0.011450             2 -117.687
## 46 -1.9610 -0.15650          0.253100 -0.7045            0.0329200  5 -110.372
## 25  4.0030                            -0.2759 -0.010390             3 -116.120
## 18  4.7460 -0.13710                           -0.012680             3 -116.207
## 49  1.0230                                    -0.011210  0.0179800  3 -117.165
## 58  8.3710 -0.29390                   -0.5758 -0.011850 -0.0139900  5 -111.985
## 37 -6.0330                   0.184000                    0.0508900  3 -117.673
## 54 -1.9280 -0.05005          0.162700         -0.011980  0.0338800  5 -112.455
## 57  1.4110                            -0.2692 -0.010350  0.0167900  4 -115.682
## 10  4.7450 -0.22870                   -0.5583                       3 -118.135
## 50  3.8810 -0.12640                           -0.012490  0.0051020  4 -116.175
## 5   2.0130                   0.141900                               2 -121.014
## 9   3.1270                            -0.3322                       2 -121.516
## 38 -6.7250  0.02906          0.189600                    0.0540500  4 -117.618
## 42  5.0310 -0.23270                   -0.5646           -0.0016620  4 -118.132
## 6   2.4260 -0.06153          0.134100                               3 -120.707
## 41 -0.4480                            -0.3019            0.0230600  3 -120.708
## 1   2.7780                                                          1 -124.240
## 33 -1.6440                                               0.0287800  2 -122.841
## 2   3.3320 -0.09176                                                 2 -123.498
## 34 -0.6381 -0.05366                                      0.0243400  3 -122.623
##     AICc  delta weight
## 28 152.5   0.00  0.525
## 27 153.1   0.56  0.397
## 31 158.6   6.08  0.025
## 59 159.1   6.56  0.020
## 32 160.1   7.58  0.012
## 60 160.8   8.28  0.008
## 19 161.4   8.89  0.006
## 20 162.4   9.90  0.004
## 23 165.3  12.75  0.001
## 51 165.3  12.81  0.001
## 63 166.3  13.81  0.001
## 24 167.9  15.43  0.000
## 52 168.7  16.19  0.000
## 55 171.3  18.75  0.000
## 11 171.6  19.04  0.000
## 64 173.2  20.74  0.000
## 12 173.4  20.89  0.000
## 43 173.9  21.43  0.000
## 15 174.9  22.35  0.000
## 56 176.5  23.95  0.000
## 47 177.2  24.70  0.000
## 16 178.1  25.55  0.000
## 44 178.9  26.36  0.000
## 35 182.5  29.96  0.000
## 3  184.8  32.26  0.000
## 48 184.9  32.42  0.000
## 36 187.1  34.55  0.000
## 39 187.2  34.65  0.000
## 7  187.9  35.42  0.000
## 4  188.3  35.74  0.000
## 8  192.4  39.93  0.000
## 40 193.3  40.79  0.000
## 30 231.2  78.65  0.000
## 29 233.2  80.74  0.000
## 61 234.4  81.94  0.000
## 14 236.1  83.60  0.000
## 45 236.7  84.15  0.000
## 13 237.4  84.94  0.000
## 26 238.1  85.59  0.000
## 21 238.1  85.60  0.000
## 53 238.9  86.38  0.000
## 62 239.1  86.63  0.000
## 22 240.7  88.18  0.000
## 17 240.7  88.20  0.000
## 46 240.7  88.23  0.000
## 25 241.2  88.73  0.000
## 18 241.4  88.91  0.000
## 49 243.3  90.82  0.000
## 58 244.0  91.46  0.000
## 37 244.3  91.84  0.000
## 54 244.9  92.40  0.000
## 57 245.1  92.57  0.000
## 10 245.3  92.76  0.000
## 50 246.1  93.55  0.000
## 5  247.4  94.85  0.000
## 9  248.4  95.86  0.000
## 38 249.0  96.44  0.000
## 42 250.0  97.47  0.000
## 6  250.4  97.90  0.000
## 41 250.4  97.91  0.000
## 1  250.9  98.37  0.000
## 33 251.0  98.51  0.000
## 2  252.3  99.82  0.000
## 34 254.2 101.74  0.000
## Models ranked by AICc(x)
summary(model.avg(model_dredge, subset = delta <= 2))
## 
## Call:
## model.avg(object = model_dredge, subset = delta <= 2)
## 
## Component model call: 
## glm(formula = NumWebs ~ <2 unique rhs>, family = poisson, data = 
##      sites_center)
## 
## Component models: 
##      df logLik   AICc delta weight
## 1234  5 -66.25 152.51  0.00   0.57
## 234   4 -69.68 153.07  0.56   0.43
## 
## Term codes: 
##     log_patch       log_tdh  log_TotalSub road_length_m 
##             1             2             3             4 
## 
## Model-averaged coefficients:  
## (full average) 
##                Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   12.284984   1.481836    1.650697   7.442  < 2e-16 ***
## log_patch     -0.167057   0.168698    0.178365   0.937 0.348965    
## log_tdh       -1.560993   0.192352    0.229029   6.816  < 2e-16 ***
## log_TotalSub  -0.659231   0.186648    0.218846   3.012 0.002593 ** 
## road_length_m -0.021230   0.004766    0.005608   3.786 0.000153 ***
##  
## (conditional average) 
##                Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   12.284984   1.481836    1.650697   7.442  < 2e-16 ***
## log_patch     -0.293327   0.113716    0.137194   2.138 0.032513 *  
## log_tdh       -1.560993   0.192352    0.229029   6.816  < 2e-16 ***
## log_TotalSub  -0.659231   0.186648    0.218846   3.012 0.002593 ** 
## road_length_m -0.021230   0.004766    0.005608   3.786 0.000153 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# We'll use the conditional average
#TotalSub
-exp(0.659231)
## [1] -1.933305
exp(0.186648)
## [1] 1.205203
#patch_area_km
-exp(0.293327)
## [1] -1.340881
exp(0.113716)
## [1] 1.120434
#Traffic_Dist_Highway
-exp(1.560993)
## [1] -4.763549
exp(0.192352)
## [1] 1.212097
#road_length_m
-0.021230 
## [1] -0.02123
0.004766
## [1] 0.004766
car::vif(get.models(model_dredge, 1)[[1]]) # looking for values < 2
##     log_patch       log_tdh  log_TotalSub road_length_m 
##      1.306159      1.088092      1.191050      1.174118

Overall Coefficients: Estimate Std. Error z value Pr(>|z|)
(Intercept) 7.55449 0.87254 8.658 < 2e-16 log_tdh -1.26833 0.19726 -6.430 1.28e-10 negative log_tdr 0.18726 0.04729 3.960 7.49e-05 *** positive log_TotalSub -1.02619 0.09768 -10.505 < 2e-16 *** negative

R-squared = 0.556197

Urban Forest Coefficients: Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.0986 0.1826 6.017 1.77e-09 ***

R-squared = -6.661338e-16

Urban Center Coefficients: Estimate Std. Error z value Pr(>|z|)
(Intercept) 13.145965 1.236110 10.635 < 2e-16 log_patch -0.293327 0.113716 -2.579 0.0099 negative log_tdh -1.539660 0.190172 -8.096 5.67e-16 ** negative log_TotalSub -0.717109 0.183835 -3.901 9.59e-05 *** negative road_length_m -0.022553 0.004713 -4.785 1.71e-06 *** negative

R-squared = 0.5814882

(conditional average) Estimate Std. Error Adjusted SE z value Pr(>|z|)
(Intercept) 12.284984 1.481836 1.650697 7.442 < 2e-16 log_patch -0.293327 0.113716 0.137194 2.138 0.032513 negative log_tdh -1.560993 0.192352 0.229029 6.816 < 2e-16 * negative log_TotalSub -0.659231 0.186648 0.218846 3.012 0.002593 ** negative road_length_m -0.021230 0.004766 0.005608 3.786 0.000153 *** negative

Number of Webs by Land

webs_land <- glm(NumWebs ~ Land, data = sites, family = poisson)
summary(webs_land)
## 
## Call:
## glm(formula = NumWebs ~ Land, family = poisson, data = sites)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -5.0080  -2.2943  -0.2411   0.2001  10.2599  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          1.1451     0.2132   5.371 7.82e-08 ***
## LandWoody Wetlands  -0.1643     0.4129  -0.398    0.691    
## LandUrban, High      1.5290     0.2504   6.107 1.02e-09 ***
## LandUrban, Medium    1.6481     0.2359   6.986 2.83e-12 ***
## LandUrban, Low       1.7726     0.2692   6.584 4.57e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 305.94  on 21  degrees of freedom
## Residual deviance: 199.27  on 17  degrees of freedom
## AIC: 288.03
## 
## Number of Fisher Scoring iterations: 6
car::Anova(webs_land, test.statistic = "LR")
## Analysis of Deviance Table (Type II tests)
## 
## Response: NumWebs
##      LR Chisq Df Pr(>Chisq)    
## Land   106.67  4  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
comp <- glht(webs_land, linfct = mcp(Land = "Tukey"))
summary(comp)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: glm(formula = NumWebs ~ Land, family = poisson, data = sites)
## 
## Linear Hypotheses:
##                                        Estimate Std. Error z value Pr(>|z|)    
## Woody Wetlands - Deciduous Forest == 0  -0.1643     0.4129  -0.398    0.994    
## Urban, High - Deciduous Forest == 0      1.5290     0.2504   6.107   <1e-04 ***
## Urban, Medium - Deciduous Forest == 0    1.6481     0.2359   6.986   <1e-04 ***
## Urban, Low - Deciduous Forest == 0       1.7726     0.2692   6.584   <1e-04 ***
## Urban, High - Woody Wetlands == 0        1.6933     0.3771   4.490   <1e-04 ***
## Urban, Medium - Woody Wetlands == 0      1.8124     0.3677   4.929   <1e-04 ***
## Urban, Low - Woody Wetlands == 0         1.9369     0.3899   4.968   <1e-04 ***
## Urban, Medium - Urban, High == 0         0.1191     0.1657   0.719    0.948    
## Urban, Low - Urban, High == 0            0.2436     0.2104   1.158    0.760    
## Urban, Low - Urban, Medium == 0          0.1246     0.1930   0.646    0.965    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
sites <- sites %>% 
  mutate(Land = fct_relevel(Land, "Deciduous Forest", 
                            "Woody Wetlands", 
                            "Urban, High", 
                            "Urban, Medium", 
                            "Urban, Low"))
nd <- data.frame(Land = factor(levels(sites$Land)))
predictions <- augment(webs_land, newdata = nd, se_fit = TRUE, type.predict = "response")
predictions <- predictions %>% 
  rename("rate" = ".fitted", "SE" = ".se.fit")
#predictions <- summary(emmeans(webs_land, ~Land),type = "response")

ggplot(sites, aes(x = Land, y = NumWebs)) + 
  geom_jitter(color = "grey", width = 0.1, size = 0.5) +
  geom_point(aes(x = Land, y = rate, color = Land), size = 1, data = predictions) + 
  geom_errorbar(aes(x = Land, 
                  ymin = rate - SE, 
                  ymax = rate + SE, color = Land), data = predictions, inherit.aes = FALSE, width = 0.25, size = 1) +
  theme_classic() +
  scale_color_manual(values = c("Deciduous Forest" = "#1b9e77", "Woody Wetlands" = "#1b9e77", "Urban, High" = "#d95f02", "Urban, Medium" = "#d95f02", "Urban, Low" = "#d95f02")) +
  xlab("Land Cover Class") +
  ylab("Number of Webs") + 
  theme(text = element_text(size = 10)) + 
  theme(axis.text.x=element_text(colour="black", size=10)) + 
  theme(axis.text.y=element_text(colour="black", size=10)) +
  theme(legend.position = "none") +
  theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust=1)) +
  annotate(geom="text", x=1, y=max(sites$NumWebs), label="A", color="black", size = 3) +
  annotate(geom="text", x=2, y=max(sites$NumWebs), label="A", color="black", size = 3) +
  annotate(geom="text", x=3, y=max(sites$NumWebs), label="B", color="black", size = 3) +
  annotate(geom="text", x=4, y=max(sites$NumWebs), label="B", color="black", size = 3) +
  annotate(geom="text", x=5, y=max(sites$NumWebs), label="B", color="black", size = 3)

The number of webs significantly varies by land cover class (chi = 106.67, df = 4, p < 0.001).

Response Variable: Number of Spiders

Number of Spiders by Location

ggarrange(ggplot(sites, aes(x = NumSpiders)) +
            geom_density(),
          ggplot(sites, aes(x = log(NumSpiders))) +
            geom_density(), ncol = 2)

shapiro.test(sites$NumSpiders) # not normal
## 
##  Shapiro-Wilk normality test
## 
## data:  sites$NumSpiders
## W = 0.66197, p-value = 6.667e-06
shapiro.test(log(sites$NumSpiders)) # not normal
## 
##  Shapiro-Wilk normality test
## 
## data:  log(sites$NumSpiders)
## W = 0.87314, p-value = 0.008961
spiders <- glm(NumSpiders ~ Location, data = sites, family = poisson)
summary(spiders)
## 
## Call:
## glm(formula = NumSpiders ~ Location, family = poisson, data = sites)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.7199  -0.9946  -0.0795   0.1385   4.7173  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            0.6419     0.2294   2.798  0.00515 ** 
## LocationUrban Center   1.2427     0.2555   4.863 1.15e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 104.356  on 21  degrees of freedom
## Residual deviance:  75.017  on 20  degrees of freedom
## AIC: 142.87
## 
## Number of Fisher Scoring iterations: 5
spiders2 <- glm(log(NumSpiders) ~ Location, data = sites)
summary(spiders2)
## 
## Call:
## glm(formula = log(NumSpiders) ~ Location, data = sites)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3213  -0.5663   0.1268   0.5169   1.7698  
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)  
## (Intercept)            0.5663     0.2870   1.973   0.0624 .
## LocationUrban Center   0.7550     0.3886   1.943   0.0662 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.8235385)
## 
##     Null deviance: 19.580  on 21  degrees of freedom
## Residual deviance: 16.471  on 20  degrees of freedom
## AIC: 62.065
## 
## Number of Fisher Scoring iterations: 2
# Both methods yield similar results
sites %>% 
  group_by(Location) %>% 
  summarize(mean = mean(NumSpiders),
            se = plotrix::std.error(NumSpiders))
## # A tibble: 2 × 3
##   Location      mean    se
##   <fct>        <dbl> <dbl>
## 1 Urban Forest  1.9  0.233
## 2 Urban Center  6.58 1.99

The number of spiders present by site is signficantly higher in the urban center (z = 4.863, df = 1, 20, p < 0.001).

Number of Spiders by Predictors

# This is one way to look at the overall data, but I think I end up not using this information
models = list(NumSpiders ~ 1,
              NumSpiders ~ log_TotalSub, 
              NumSpiders ~ log_tdr, 
              NumSpiders ~ log_tdh,
              NumSpiders ~ log_TotalSub + log_tdr,
              NumSpiders ~ log_TotalSub + log_tdh,
              NumSpiders ~ log_tdr + log_tdh,
              NumSpiders ~ log_TotalSub + log_tdr + log_tdh)
fits = lapply(models, glm, data = sites, family = "poisson")
modnames = sapply(models, function(ff)deparse(ff[[3]]))
pander(aictab(fits, modname = modnames), caption="Table 1. AICc model selection table for Search Distance to the Focal Web Overall.", split.tables = Inf)
Table 1. AICc model selection table for Search Distance to the Focal Web Overall.
  Modnames K AICc Delta_AICc ModelLik AICcWt LL Cum.Wt
6 log_TotalSub + log_tdh 3 119.6 0 1 0.8187 -56.12 0.8187
8 log_TotalSub + log_tdr + log_tdh 4 122.6 3.015 0.2215 0.1813 -56.12 1
2 log_TotalSub 2 139.9 20.3 3.898e-05 3.191e-05 -67.62 1
5 log_TotalSub + log_tdr 3 142.1 22.52 1.288e-05 1.055e-05 -67.38 1
4 log_tdh 2 168.4 48.82 2.507e-11 2.052e-11 -81.88 1
7 log_tdr + log_tdh 3 169.5 49.92 1.445e-11 1.183e-11 -81.08 1
1 1 1 170.4 50.84 9.146e-12 7.487e-12 -84.1 1
3 log_tdr 2 170.5 50.92 8.779e-12 7.187e-12 -82.93 1
# The top model keeps all 3 variables and no other model is within 2 Model Likelihood Units
summary(glm(NumSpiders ~ log_TotalSub + log_tdr + log_tdh, data = sites, family = poisson))
## 
## Call:
## glm(formula = NumSpiders ~ log_TotalSub + log_tdr + log_tdh, 
##     family = poisson, data = sites)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5806  -1.2639  -0.3472   0.5572   3.2858  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   7.360057   1.229989   5.984 2.18e-09 ***
## log_TotalSub -0.953104   0.137313  -6.941 3.89e-12 ***
## log_tdr      -0.004858   0.071234  -0.068    0.946    
## log_tdh      -1.196433   0.282692  -4.232 2.31e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 104.356  on 21  degrees of freedom
## Residual deviance:  48.382  on 18  degrees of freedom
## AIC: 120.23
## 
## Number of Fisher Scoring iterations: 5
global <- glm(NumSpiders ~ log_TotalSub + log_tdr + log_tdh, data = sites, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_ns_o <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_ns_o)
out
## 
## Call:
## glm(formula = NumSpiders ~ log_tdh + log_TotalSub + 1, family = poisson, 
##     data = sites)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5733  -1.2492  -0.3432   0.5522   3.2489  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    7.3234     1.1068   6.617 3.67e-11 ***
## log_tdh       -1.1932     0.2785  -4.284 1.84e-05 ***
## log_TotalSub  -0.9523     0.1370  -6.952 3.61e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 104.356  on 21  degrees of freedom
## Residual deviance:  48.387  on 19  degrees of freedom
## AIC: 118.24
## 
## Number of Fisher Scoring iterations: 5
info <- out[[12]]
#TotalSub (log_TotalSub)
exp(info[3, 1])
## [1] 0.3858395
exp(info[3, 2])
## [1] 1.146825
#Traffic_Dist_Highway (log_tdh)
exp(info[2, 1])
## [1] 0.3032556
exp(info[2, 2])
## [1] 1.321175
with(summary(top_model_ns_o), 1 - deviance/null.deviance) # This gives the R squared value
## [1] 0.5363275
#summary(model.avg(model_dredge, subset = delta <= 2))
#summary(model.avg(model_dredge))
model.sel(model_dredge) #estimates same signs
## Global model call: glm(formula = NumSpiders ~ log_TotalSub + log_tdr + log_tdh, 
##     family = poisson, data = sites)
## ---
## Model selection table 
##    (Int) log_tdh   log_tdr log_TtS df  logLik  AICc delta weight
## 6 7.3230 -1.1930           -0.9523  3 -56.119 119.6  0.00  0.819
## 8 7.3600 -1.1960 -0.004858 -0.9531  4 -56.116 122.6  3.01  0.181
## 5 2.5440                   -0.7056  2 -67.622 139.9 20.30  0.000
## 7 2.2760          0.053410 -0.6989  3 -67.378 142.1 22.52  0.000
## 2 3.0700 -0.4256                    2 -81.879 168.4 48.82  0.000
## 4 2.4780 -0.3834  0.091100          3 -81.079 169.5 49.92  0.000
## 1 1.4940                            1 -84.103 170.4 50.84  0.000
## 3 0.9569          0.112200          2 -82.928 170.5 50.92  0.000
## Models ranked by AICc(x)
car::vif(get.models(model_dredge, 1)[[1]]) #looking for less than 2
##      log_tdh log_TotalSub 
##     1.237425     1.237425
# urban forest
global <- glm(NumSpiders ~ tree100m + log_TotalSub + log_tdh + log_patch + log_tdr, data = sites_forest, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_ns_f <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_ns_f)
out
## 
## Call:
## glm(formula = NumSpiders ~ 1, family = poisson, data = sites_forest)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.71853  -0.52092   0.07192   0.07192   0.73522  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)   0.6419     0.2294   2.798  0.00515 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2.6558  on 9  degrees of freedom
## Residual deviance: 2.6558  on 9  degrees of freedom
## AIC: 29.708
## 
## Number of Fisher Scoring iterations: 4
info <- out[[12]]
with(summary(top_model_ns_f), 1 - deviance/null.deviance) # This gives the R squared value
## [1] -8.881784e-16
model.sel(model_dredge) #estimates same signs
## Global model call: glm(formula = NumSpiders ~ tree100m + log_TotalSub + log_tdh + 
##     log_patch + log_tdr, family = poisson, data = sites_forest)
## ---
## Model selection table 
##     (Int)  log_ptc  log_tdh log_tdr log_TtS      t10 df  logLik AICc delta
## 1  0.6419                                             1 -13.854 30.2  0.00
## 5  1.1250                   -0.1238                   2 -13.708 33.1  2.92
## 17 1.1650                                   -0.01122  2 -13.742 33.2  2.99
## 2  1.7930 -0.10970                                    2 -13.802 33.3  3.11
## 9  0.5202                           0.04853           2 -13.850 33.4  3.21
## 3  0.7653          -0.03509                           2 -13.853 33.4  3.21
## 21 2.8130                   -0.2448         -0.02608  3 -13.312 36.6  6.42
## 6  2.5900 -0.13370          -0.1395                   3 -13.630 37.3  7.05
## 13 0.7519                   -0.1540 0.19600           3 -13.654 37.3  7.10
## 7  2.0200          -0.22290 -0.1522                   3 -13.665 37.3  7.12
## 18 2.1260 -0.09493                          -0.01046  3 -13.705 37.4  7.20
## 19 1.2380          -0.02173                 -0.01114  3 -13.742 37.5  7.28
## 25 1.1190                           0.01673 -0.01113  3 -13.742 37.5  7.28
## 4  2.0800 -0.11450 -0.06705                           3 -13.797 37.6  7.39
## 10 1.6710 -0.11070                  0.05309           3 -13.798 37.6  7.39
## 11 0.6184          -0.02513         0.04462           3 -13.849 37.7  7.49
## 23 4.7680          -0.43260 -0.3210         -0.02898  4 -13.167 42.3 12.13
## 29 2.5070                   -0.3123 0.30140 -0.03008  4 -13.183 42.4 12.16
## 22 3.5930 -0.08358          -0.2420         -0.02421  4 -13.284 42.6 12.36
## 8  4.0810 -0.15930 -0.30060 -0.1815                   4 -13.556 43.1 12.90
## 14 2.2610 -0.14090          -0.1716 0.21130           4 -13.568 43.1 12.93
## 15 1.7070          -0.24150 -0.1892 0.20860           4 -13.605 43.2 13.00
## 20 2.3870 -0.10100 -0.05849                 -0.01028  4 -13.701 43.4 13.19
## 26 2.0640 -0.09647                  0.02886 -0.01033  4 -13.704 43.4 13.20
## 27 1.1890          -0.01807         0.01317 -0.01108  4 -13.741 43.5 13.27
## 12 1.9360 -0.11470 -0.05711         0.04423           4 -13.794 43.6 13.38
## 31 4.6300          -0.47640 -0.4083 0.32910 -0.03313  5 -13.017 51.0 20.83
## 24 6.4050 -0.13430 -0.51230 -0.3320         -0.02689  5 -13.096 51.2 20.98
## 30 3.3740 -0.09490          -0.3101 0.31320 -0.02812  5 -13.147 51.3 21.09
## 16 3.8620 -0.16930 -0.32710 -0.2241 0.23250           5 -13.484 52.0 21.76
## 28 2.3240 -0.10150 -0.05338         0.01888 -0.01021  5 -13.701 52.4 22.19
## 32 6.5900 -0.15850 -0.58480 -0.4312 0.36300 -0.03117  6 -12.922 65.8 35.64
##    weight
## 1   0.423
## 5   0.098
## 17  0.095
## 2   0.089
## 9   0.085
## 3   0.085
## 21  0.017
## 6   0.012
## 13  0.012
## 7   0.012
## 18  0.012
## 19  0.011
## 25  0.011
## 4   0.011
## 10  0.011
## 11  0.010
## 23  0.001
## 29  0.001
## 22  0.001
## 8   0.001
## 14  0.001
## 15  0.001
## 20  0.001
## 26  0.001
## 27  0.001
## 12  0.001
## 31  0.000
## 24  0.000
## 30  0.000
## 16  0.000
## 28  0.000
## 32  0.000
## Models ranked by AICc(x)
#summary(model.avg(model_dredge, subset = delta <= 2))
#car::vif(get.models(model_dredge, 1)[[1]]) #looking for less than 2

# urban center
global <- glm(NumSpiders ~ log_TotalSub + log_tdr + log_tdh + spec_rad + log_patch + road_length_m, data = sites_center, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_ns_c <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_ns_c)
out
## 
## Call:
## glm(formula = NumSpiders ~ log_patch + log_tdh + log_TotalSub + 
##     road_length_m + 1, family = poisson, data = sites_center)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.81882  -0.84088  -0.06627   0.61185   1.84897  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   14.516312   2.195703   6.611 3.81e-11 ***
## log_patch     -0.618055   0.205145  -3.013 0.002589 ** 
## log_tdh       -1.393026   0.286317  -4.865 1.14e-06 ***
## log_TotalSub  -1.331282   0.357971  -3.719 0.000200 ***
## road_length_m -0.025605   0.007753  -3.303 0.000958 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 72.361  on 11  degrees of freedom
## Residual deviance: 17.493  on  7  degrees of freedom
## AIC: 66.291
## 
## Number of Fisher Scoring iterations: 5
info <- out[[12]]
#TotalSub
exp(info[4, 1])
## [1] 0.2641384
exp(info[4, 2])
## [1] 1.430424
#patch_area_km
exp(info[2, 1])
## [1] 0.5389918
exp(info[2, 2])
## [1] 1.227703
#Traffic_Dist_Highway
exp(info[3, 1])
## [1] 0.2483228
exp(info[3, 2])
## [1] 1.331515
#road_length
info[5, 1]
## [1] -0.02560466
info[5, 2]
## [1] 0.007752731
with(summary(top_model_ns_c), 1 - deviance/null.deviance) # This gives the R squared value
## [1] 0.7582528
model.sel(model_dredge) #estimates same signs
## Global model call: glm(formula = NumSpiders ~ log_TotalSub + log_tdr + log_tdh + 
##     spec_rad + log_patch + road_length_m, family = poisson, data = sites_center)
## ---
## Model selection table 
##      (Int)  log_ptc log_tdh   log_tdr log_TtS rod_lng_m   spc_rad df  logLik
## 28 14.5200 -0.61810  -1.393           -1.3310 -0.025600            5 -28.146
## 27 10.1100           -1.486           -0.8468 -0.019220            4 -33.392
## 12 10.5300 -0.41940  -1.305           -1.1540                      4 -34.818
## 59  1.3070           -1.551           -0.8344 -0.019100  0.058520  5 -31.812
## 32 14.5900 -0.59230  -1.420 -0.054740 -1.2230 -0.024950            6 -27.922
## 60 12.5000 -0.58610  -1.415           -1.3040 -0.025240  0.011970  6 -28.093
## 43 -3.0470           -1.575           -0.7302            0.075560  4 -35.661
## 23 10.1500           -1.443 -0.194200         -0.019590            4 -35.785
## 11  8.1870           -1.447           -0.7967                      3 -38.200
## 31 10.6100           -1.526 -0.096830 -0.6999 -0.019050            5 -32.735
## 19  8.3150           -1.227                   -0.019570            3 -39.096
## 20 10.8900 -0.34340  -1.230                   -0.025330            4 -36.907
## 24 12.3200 -0.32160  -1.416 -0.185900         -0.023900            5 -33.803
## 51 -0.7455           -1.298                   -0.017780  0.059480  4 -37.174
## 44  3.0050 -0.32240  -1.429           -1.0260            0.047220  5 -34.041
## 16 10.9900 -0.40750  -1.355 -0.086330 -1.0410                      5 -34.340
## 15  8.8330           -1.509 -0.103000 -0.6735                      4 -37.537
## 35 -5.7790           -1.325                              0.082410  3 -40.533
## 55  5.6420           -1.441 -0.160300         -0.018520  0.027430  5 -35.440
## 47 -2.1070           -1.582 -0.027620 -0.6988            0.070380  5 -35.621
## 7   8.3550           -1.414 -0.210400                              3 -41.193
## 39 -0.9622           -1.458 -0.138200                    0.059100  4 -39.306
## 63  2.3510           -1.555 -0.027550 -0.7928 -0.018920  0.052420  6 -31.773
## 52  4.3960 -0.24320  -1.272                   -0.022360  0.037580  5 -36.313
## 3   6.2940           -1.160                                        2 -44.518
## 8   9.1680 -0.16750  -1.354 -0.216600                              4 -40.429
## 36 -5.3840 -0.02458  -1.315                              0.080560  4 -40.519
## 56 14.8900 -0.36160  -1.414 -0.203100         -0.025020 -0.013930  6 -33.743
## 48  4.7690 -0.33230  -1.432 -0.046520 -0.9851            0.037580  6 -33.930
## 4   6.9950 -0.14440  -1.115                                        3 -43.934
## 64 15.0600 -0.59830  -1.416 -0.058250 -1.2230 -0.025010 -0.002777  7 -27.920
## 40  0.4998 -0.07391  -1.430 -0.147800                    0.052130  5 -39.183
## 26  7.6110 -0.54850                   -1.2270 -0.013420            4 -44.356
## 10  5.8300 -0.47000                   -1.1100                      3 -47.214
## 42  2.8150 -0.42810                   -1.0460            0.017520  4 -47.085
## 14  5.7870 -0.47030          0.010440 -1.1190                      4 -47.208
## 30  7.5550 -0.55460          0.028000 -1.2610 -0.013640            5 -44.308
## 58  6.5570 -0.53210                   -1.2070 -0.013350  0.006031  5 -44.340
## 25  3.3350                            -0.5552 -0.009836            3 -50.162
## 41 -6.4160                            -0.5222            0.057370  3 -50.183
## 9   2.4650                            -0.5681                      2 -52.208
## 57 -4.8090                            -0.5550 -0.010160  0.053040  4 -48.362
## 49 -5.2100                                    -0.010660  0.052310  3 -51.000
## 33 -7.7040                                               0.062290  2 -52.842
## 17  2.9170                                    -0.011470            2 -52.889
## 18  4.4140 -0.21780                           -0.013590            3 -51.369
## 45 -9.0390                   0.096520 -0.6097            0.071700  4 -49.714
## 13  2.4010                   0.015170 -0.5828                      3 -52.196
## 1   1.8850                                                         1 -55.580
## 46  1.2770 -0.41160          0.042670 -1.0600            0.025520  5 -47.003
## 29  3.2980                   0.008705 -0.5644 -0.009830            4 -50.158
## 21  3.2590                  -0.069050         -0.011270            3 -52.549
## 34 -6.3200 -0.07443                                      0.056200  3 -52.674
## 2   2.8610 -0.16250                                                2 -54.621
## 37 -7.8320                   0.004515                    0.062970  3 -52.841
## 50 -2.2580 -0.13540                           -0.011980  0.039180  4 -50.537
## 22  4.8380 -0.22440         -0.081340         -0.013130            4 -50.894
## 61 -7.7730                   0.101700 -0.6649 -0.010650  0.069830  5 -47.852
## 53 -5.1540                  -0.001917         -0.010660  0.052010  4 -50.999
## 5   2.2750                  -0.074990                              2 -55.201
## 6   3.4850 -0.18230         -0.096930                              3 -54.001
## 62  4.7150 -0.51690          0.049730 -1.2360 -0.013650  0.016070  6 -44.228
## 38 -5.7770 -0.08121         -0.014750                    0.053450  4 -52.662
## 54 -0.9145 -0.15180         -0.034010         -0.012040  0.032290  5 -50.475
##     AICc delta weight
## 28  76.3  0.00  0.784
## 27  80.5  4.21  0.096
## 12  83.3  7.06  0.023
## 59  83.6  7.33  0.020
## 32  84.6  8.35  0.012
## 60  85.0  8.69  0.010
## 43  85.0  8.74  0.010
## 23  85.3  8.99  0.009
## 11  85.4  9.11  0.008
## 31  85.5  9.18  0.008
## 19  87.2 10.90  0.003
## 20  87.5 11.24  0.003
## 24  87.6 11.31  0.003
## 51  88.1 11.77  0.002
## 44  88.1 11.79  0.002
## 16  88.7 12.39  0.002
## 15  88.8 12.50  0.002
## 35  90.1 13.77  0.001
## 55  90.9 14.59  0.001
## 47  91.2 14.95  0.000
## 7   91.4 15.09  0.000
## 39  92.3 16.03  0.000
## 63  92.3 16.06  0.000
## 52  92.6 16.33  0.000
## 3   94.4 18.08  0.000
## 8   94.6 18.28  0.000
## 36  94.8 18.46  0.000
## 56  96.3 19.99  0.000
## 48  96.7 20.37  0.000
## 4   96.9 20.58  0.000
## 64  97.8 21.55  0.000
## 40  98.4 22.07  0.000
## 26 102.4 26.13  0.000
## 10 103.4 27.14  0.000
## 42 107.9 31.59  0.000
## 14 108.1 31.84  0.000
## 30 108.6 32.32  0.000
## 58 108.7 32.39  0.000
## 25 109.3 33.03  0.000
## 41 109.4 33.08  0.000
## 9  109.7 33.46  0.000
## 57 110.4 34.15  0.000
## 49 111.0 34.71  0.000
## 33 111.0 34.73  0.000
## 17 111.1 34.82  0.000
## 18 111.7 35.45  0.000
## 45 113.1 36.85  0.000
## 13 113.4 37.10  0.000
## 1  113.6 37.27  0.000
## 46 114.0 37.72  0.000
## 29 114.0 37.74  0.000
## 21 114.1 37.81  0.000
## 34 114.3 38.06  0.000
## 2  114.6 38.28  0.000
## 37 114.7 38.39  0.000
## 50 114.8 38.50  0.000
## 22 115.5 39.21  0.000
## 61 115.7 39.41  0.000
## 53 115.7 39.42  0.000
## 5  115.7 39.44  0.000
## 6  117.0 40.71  0.000
## 62 117.3 40.97  0.000
## 38 119.0 42.75  0.000
## 54 120.9 44.66  0.000
## Models ranked by AICc(x)
#summary(model.avg(model_dredge, subset = delta <= 2))
car::vif(get.models(model_dredge, 1)[[1]]) # looking for values < 2
##     log_patch       log_tdh  log_TotalSub road_length_m 
##      1.827097      1.084576      1.573369      1.285511

Overall Coefficients: Estimate Std. Error z value Pr(>|z|)
(Intercept) 7.3234 1.1068 6.617 3.67e-11 log_tdh -1.1932 0.2785 -4.284 1.84e-05 negative log_TotalSub -0.9523 0.1370 -6.952 3.61e-12 *** negative

R-squared = 0.5363275

Urban Forest Coefficients: Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.6419 0.2294 2.798 0.00515 **

R-squared = -8.881784e-16

Urban Center Coefficients: Estimate Std. Error z value Pr(>|z|)
(Intercept) 14.516312 2.195703 6.611 3.81e-11 log_patch -0.618055 0.205145 -3.013 0.002589 negative log_tdh -1.393026 0.286317 -4.865 1.14e-06 ** negative log_TotalSub -1.331282 0.357971 -3.719 0.000200 *** negative road_length_m -0.025605 0.007753 -3.303 0.000958 *** negative

R-squared = 0.7582528

Number of Spiders by Land

spiders_land <- glm(NumSpiders ~ Land, data = sites, family = poisson)
summary(spiders_land)
## 
## Call:
## glm(formula = NumSpiders ~ Land, family = poisson, data = sites)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.9950  -0.8393  -0.2352   0.3782   5.1978  
## 
## Coefficients:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          0.7621     0.2582   2.952 0.003160 ** 
## LandWoody Wetlands  -0.4745     0.5627  -0.843 0.399154    
## LandUrban, High      1.2528     0.3162   3.962 7.45e-05 ***
## LandUrban, Medium    0.9725     0.3100   3.137 0.001705 ** 
## LandUrban, Low       1.2528     0.3651   3.431 0.000602 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 104.356  on 21  degrees of freedom
## Residual deviance:  72.705  on 17  degrees of freedom
## AIC: 146.56
## 
## Number of Fisher Scoring iterations: 5
car::Anova(spiders_land, test.statistic = "LR")
## Analysis of Deviance Table (Type II tests)
## 
## Response: NumSpiders
##      LR Chisq Df Pr(>Chisq)    
## Land    31.65  4  2.255e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
comp <- glht(spiders_land, linfct = mcp(Land = "Tukey"))
summary(comp)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: glm(formula = NumSpiders ~ Land, family = poisson, data = sites)
## 
## Linear Hypotheses:
##                                          Estimate Std. Error z value Pr(>|z|)
## Woody Wetlands - Deciduous Forest == 0 -4.745e-01  5.627e-01  -0.843  0.91161
## Urban, High - Deciduous Forest == 0     1.253e+00  3.162e-01   3.962  < 0.001
## Urban, Medium - Deciduous Forest == 0   9.725e-01  3.100e-01   3.137  0.01339
## Urban, Low - Deciduous Forest == 0      1.253e+00  3.651e-01   3.431  0.00482
## Urban, High - Woody Wetlands == 0       1.727e+00  5.323e-01   3.245  0.00907
## Urban, Medium - Woody Wetlands == 0     1.447e+00  5.286e-01   2.737  0.04467
## Urban, Low - Woody Wetlands == 0        1.727e+00  5.627e-01   3.069  0.01644
## Urban, Medium - Urban, High == 0       -2.803e-01  2.505e-01  -1.119  0.78534
## Urban, Low - Urban, High == 0           2.820e-14  3.162e-01   0.000  1.00000
## Urban, Low - Urban, Medium == 0         2.803e-01  3.100e-01   0.904  0.88867
##                                           
## Woody Wetlands - Deciduous Forest == 0    
## Urban, High - Deciduous Forest == 0    ***
## Urban, Medium - Deciduous Forest == 0  *  
## Urban, Low - Deciduous Forest == 0     ** 
## Urban, High - Woody Wetlands == 0      ** 
## Urban, Medium - Woody Wetlands == 0    *  
## Urban, Low - Woody Wetlands == 0       *  
## Urban, Medium - Urban, High == 0          
## Urban, Low - Urban, High == 0             
## Urban, Low - Urban, Medium == 0           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
sites <- sites %>% 
  mutate(Land = fct_relevel(Land, "Deciduous Forest", 
                            "Woody Wetlands", 
                            "Urban, High", 
                            "Urban, Medium", 
                            "Urban, Low"))
nd <- data.frame(Land = factor(levels(sites$Land)))
predictions <- augment(spiders_land, newdata = nd, se_fit = TRUE, type.predict = "response")
predictions <- predictions %>% 
  rename("rate" = ".fitted", "SE" = ".se.fit")
#predictions <- summary(emmeans(spiders_land, ~Land),type = "response")

ggplot(sites, aes(x = Land, y = NumSpiders)) + 
  geom_jitter(color = "grey", width = 0.1, size = 0.5) +
  geom_point(aes(x = Land, y = rate, color = Land), size = 1, data = predictions) + 
  geom_errorbar(aes(x = Land, 
                  ymin = rate - SE, 
                  ymax = rate + SE, color = Land), data = predictions, inherit.aes = FALSE, width = 0.25, size = 1) +
  theme_classic() +
  scale_color_manual(values = c("Deciduous Forest" = "#1b9e77", "Woody Wetlands" = "#1b9e77", "Urban, High" = "#d95f02", "Urban, Medium" = "#d95f02", "Urban, Low" = "#d95f02")) +
  xlab("Land Cover Class") +
  ylab("Number of Spiders") + 
  theme(text = element_text(size = 10)) + 
  theme(axis.text.x=element_text(colour="black", size=10)) + 
  theme(axis.text.y=element_text(colour="black", size=10)) +
  theme(legend.position = "none") +
  theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust=1)) +
  annotate(geom="text", x=1, y=max(sites$NumSpiders), label="A", color="black", size = 3) +
  annotate(geom="text", x=2, y=max(sites$NumSpiders), label="A", color="black", size = 3) +
  annotate(geom="text", x=3, y=max(sites$NumSpiders), label="B", color="black", size = 3) +
  annotate(geom="text", x=4, y=max(sites$NumSpiders), label="B", color="black", size = 3) +
  annotate(geom="text", x=5, y=max(sites$NumSpiders), label="B", color="black", size = 3)

The number of spiders per site significantly varies across land cover class (Chi = 31.65, df = 4, p < 0.001).

Response Variable: Web Distance

Web Distance by Location

ggarrange(ggplot(web_dist, aes(x = RetreatDist)) +
            geom_density(),
          ggplot(web_dist, aes(x = log(RetreatDist))) +
            geom_density(), ncol = 2)

shapiro.test(web_dist$RetreatDist) # not normal
## 
##  Shapiro-Wilk normality test
## 
## data:  web_dist$RetreatDist
## W = 0.8652, p-value = 0.0002559
shapiro.test(log(web_dist$RetreatDist)) # not normal
## 
##  Shapiro-Wilk normality test
## 
## data:  log(web_dist$RetreatDist)
## W = 0.92787, p-value = 0.01536
dist <- glmer(round(RetreatDist, 0) ~ Location * factor(Neighbor) + (1 | ID), data = web_dist, family = poisson)
summary(dist)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: round(RetreatDist, 0) ~ Location * factor(Neighbor) + (1 | ID)
##    Data: web_dist
## 
##      AIC      BIC   logLik deviance df.resid 
##   1479.0   1487.3   -734.5   1469.0       34 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -15.645  -1.728  -0.013   2.028  11.418 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 0.5629   0.7503  
## Number of obs: 39, groups:  ID, 21
## 
## Fixed effects:
##                                        Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                             6.09872    0.23780  25.647  < 2e-16 ***
## LocationUrban Center                   -1.84667    0.32996  -5.597 2.18e-08 ***
## factor(Neighbor)2                       0.33655    0.02297  14.654  < 2e-16 ***
## LocationUrban Center:factor(Neighbor)2  0.41493    0.04208   9.860  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) LctnUC fc(N)2
## LctnUrbnCnt -0.721              
## fctr(Nghb)2 -0.039  0.028       
## LctnUC:(N)2  0.022 -0.076 -0.546
dist2 <- lmer(log(RetreatDist) ~ Location * factor(Neighbor) + (1 | ID), data = web_dist)
coefs <- data.frame(coef(summary(dist2)))
coefs$p.z <- 2 * (1 - pnorm(abs(coefs$t.value)))
coefs
##                                         Estimate Std..Error    t.value
## (Intercept)                             6.088014  0.2721524 22.3698674
## LocationUrban Center                   -1.998638  0.3760326 -5.3150650
## factor(Neighbor)2                       0.329375  0.3580922  0.9198049
## LocationUrban Center:factor(Neighbor)2  0.485294  0.4679854  1.0369853
##                                                 p.z
## (Intercept)                            0.000000e+00
## LocationUrban Center                   1.066193e-07
## factor(Neighbor)2                      3.576747e-01
## LocationUrban Center:factor(Neighbor)2 2.997427e-01
# Both methods yield similar results
web_dist %>% 
  group_by(Location) %>% 
  summarize(mean = mean(RetreatDist),
            se = plotrix::std.error(RetreatDist))
## # A tibble: 2 × 3
##   Location      mean    se
##   <fct>        <dbl> <dbl>
## 1 Urban Forest  568.  63.0
## 2 Urban Center  168.  48.1

Spider neighbors were significantly closer in the urban center than the urban forest (z = -5.597, df = 1, 21, p < 0.001). The second neighbor was always further from the focal web than the second (z = 14.654, df = 1, 21, p < 0.001), as this is how we defined the first and second neighbor. We also found a significant interaction between location and neighbor (z = 9.860, df = 1, 21, p < 0.001) - the distance between the first and second neighbor differed to a higher degree in the urban center than the urban forest.

Web Distance by Predictors

To look at the predictors, we will restrict to just the nearest neighbor instead of the nearest and second nearest neighbor.

#overall
models = list(round(RetreatDist,0) ~ 1,
              round(RetreatDist,0) ~ log_TotalSub,
              round(RetreatDist,0) ~ log_tdr, 
              round(RetreatDist,0) ~ log_tdh,
              round(RetreatDist,0) ~ log_TotalSub + log_tdr,
              round(RetreatDist,0) ~ log_TotalSub + log_tdh,
              round(RetreatDist,0) ~ log_tdr + log_tdh,
              round(RetreatDist,0) ~ log_TotalSub + log_tdr + log_tdh)
fits = lapply(models, glm, data = near, family = "poisson")
modnames = sapply(models, function(ff)deparse(ff[[3]]))
pander(aictab(fits, modname = modnames), caption="Table 1. AICc model selection table for Search Distance to the Focal Web.", split.tables = Inf)
Table 1. AICc model selection table for Search Distance to the Focal Web.
  Modnames K AICc Delta_AICc ModelLik AICcWt LL Cum.Wt
5 log_TotalSub + log_tdr 3 4029 0 1 0.7946 -2011 0.7946
8 log_TotalSub + log_tdr + log_tdh 4 4031 2.706 0.2584 0.2054 -2010 1
6 log_TotalSub + log_tdh 3 4353 324.5 3.488e-71 2.771e-71 -2173 1
2 log_TotalSub 2 4364 335.3 1.567e-73 1.245e-73 -2180 1
7 log_tdr + log_tdh 3 5082 1053 1.811e-229 1.439e-229 -2537 1
3 log_tdr 2 5333 1304 6.089e-284 4.838e-284 -2664 1
4 log_tdh 2 5839 1811 0 0 -2917 1
1 1 1 6042 2013 0 0 -3020 1
global <- glm(round(RetreatDist, 0) ~ log_TotalSub + log_tdr + log_tdh, data = near, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_wd_o <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_wd_o)
out
## 
## Call:
## glm(formula = round(RetreatDist, 0) ~ log_tdr + log_TotalSub + 
##     1, family = poisson, data = near)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -17.736  -12.320   -4.742    4.158   29.372  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   5.38311    0.06200   86.83   <2e-16 ***
## log_tdr      -0.18279    0.01024  -17.84   <2e-16 ***
## log_TotalSub  0.57484    0.01653   34.78   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 5895.3  on 20  degrees of freedom
## Residual deviance: 3876.6  on 18  degrees of freedom
## AIC: 4027.1
## 
## Number of Fisher Scoring iterations: 5
info <- out[[12]]
# TotalSub
exp(info[3, 1])
## [1] 1.776844
exp(info[3, 2])
## [1] 1.016665
#Traffic_Dist_Road
exp(info[2, 1])
## [1] 0.832942
exp(info[2, 1])
## [1] 0.832942
with(summary(top_model_wd_o), 1 - deviance/null.deviance)
## [1] 0.3424213
model.sel(model_dredge) #estimates same signs
## Global model call: glm(formula = round(RetreatDist, 0) ~ log_TotalSub + log_tdr + 
##     log_tdh, family = poisson, data = near)
## ---
## Model selection table 
##   (Int)  log_tdh log_tdr log_TtS df    logLik   AICc   delta weight
## 7 5.383          -0.1828  0.5748  3 -2010.560 4028.5    0.00  0.795
## 8 5.464 -0.01786 -0.1840  0.5699  4 -2010.369 4031.2    2.71  0.205
## 6 4.049  0.09758          0.6593  3 -2172.794 4353.0  324.47  0.000
## 5 4.468                   0.6308  2 -2179.572 4363.8  335.28  0.000
## 4 8.427 -0.42250 -0.2622          3 -2537.258 5081.9 1053.40  0.000
## 3 6.855          -0.2613          2 -2664.060 5332.8 1304.26  0.000
## 2 6.994 -0.34820                  2 -2917.187 5839.0 1810.51  0.000
## 1 5.701                           1 -3019.903 6042.0 2013.49  0.000
## Models ranked by AICc(x)
#summary(model.avg(model_dredge, subset = delta <= 2))
car::vif(get.models(model_dredge, 5)[[1]]) # looking for less than 2
##  log_tdh  log_tdr 
## 1.003816 1.003816
# Urban Forest
global <- glm(round(RetreatDist, 0) ~ tree100m + log_TotalSub + log_tdr + log_tdh + log_patch, data = near_forest, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_wd_f <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_wd_f)
out
## 
## Call:
## glm(formula = round(RetreatDist, 0) ~ log_tdh + log_tdr + log_TotalSub + 
##     tree100m + 1, family = poisson, data = near_forest)
## 
## Deviance Residuals: 
##       1        2        3        4        5        6        7        8  
## -10.062  -22.098   10.875   -3.513   15.652   14.285    1.742   -7.525  
##       9       10  
##  -3.027   -5.783  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   5.247233   0.253935  20.664  < 2e-16 ***
## log_tdh       0.192623   0.045552   4.229 2.35e-05 ***
## log_tdr       0.097475   0.018717   5.208 1.91e-07 ***
## log_TotalSub -0.182147   0.037627  -4.841 1.29e-06 ***
## tree100m      0.007928   0.001747   4.540 5.64e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1320.6  on 9  degrees of freedom
## Residual deviance: 1271.5  on 5  degrees of freedom
## AIC: 1360.7
## 
## Number of Fisher Scoring iterations: 5
info <- out[[12]]
#tree100m
info[5, 1]
## [1] 0.007928466
info[5, 2]
## [1] 0.001746549
#TotalSub
exp(info[4, 1])
## [1] 0.8334784
exp(info[4, 2])
## [1] 1.038344
#Road
exp(info[3, 1])
## [1] 1.102384
exp(info[3, 2])
## [1] 1.018894
#Highway
exp(info[2, 1])
## [1] 1.212425
exp(info[2, 2])
## [1] 1.046605
with(summary(top_model_wd_f), 1 - deviance/null.deviance)
## [1] 0.03720731
model.sel(model_dredge) #estimates same signs
## Global model call: glm(formula = round(RetreatDist, 0) ~ tree100m + log_TotalSub + 
##     log_tdr + log_tdh + log_patch, family = poisson, data = near_forest)
## ---
## Model selection table 
##    (Int)  log_ptc log_tdh  log_tdr log_TtS      t10 df   logLik   AICc delta
## 31 5.247          0.19260 0.097480 -0.1821 0.007928  5 -675.375 1375.7  0.00
## 29 6.163                  0.064510 -0.1803 0.005608  4 -684.326 1384.7  8.90
## 15 5.967          0.13170 0.049800 -0.1584           4 -685.602 1387.2 11.45
## 32 4.923 0.026220 0.20270 0.101200 -0.1834 0.007949  6 -674.627 1389.3 13.51
## 23 4.986          0.19020 0.063390         0.006888  4 -687.159 1390.3 14.57
## 13 6.502                  0.034750 -0.1625           3 -690.175 1390.4 14.60
## 11 6.207          0.08834          -0.1139           3 -690.668 1391.3 15.59
## 9  6.548                           -0.1257           2 -692.895 1391.5 15.76
## 25 6.413                           -0.1200 0.002561  3 -691.307 1392.6 16.86
## 30 6.071 0.008685         0.065180 -0.1807 0.005580  5 -684.239 1393.5 17.73
## 27 6.030          0.09568          -0.1078 0.002873  4 -688.740 1393.5 17.73
## 16 5.624 0.027680 0.14360 0.053540 -0.1604           5 -684.782 1394.6 18.82
## 10 6.481 0.006370                  -0.1256           3 -692.848 1395.7 19.95
## 14 6.365 0.012810         0.035870 -0.1636           4 -689.988 1396.0 20.23
## 12 6.048 0.013550 0.09262          -0.1132           4 -690.460 1396.9 21.17
## 19 5.670          0.11570                  0.003324  3 -693.669 1397.3 21.59
## 24 4.705 0.022510 0.19910 0.066360         0.006917  5 -686.614 1398.2 22.48
## 3  5.844          0.11080                            2 -696.266 1398.2 22.50
## 26 6.398 0.001467                  -0.1200 0.002548  4 -691.304 1398.6 22.86
## 7  5.649          0.13840 0.024610                   3 -694.822 1399.6 23.90
## 17 6.089                                   0.003094  2 -697.597 1400.9 25.16
## 21 5.897                  0.029460         0.004670  3 -695.805 1401.6 25.86
## 4  5.664 0.015570 0.11540                            3 -695.996 1402.0 26.24
## 28 5.940 0.008032 0.09786          -0.1074 0.002800  5 -688.668 1402.3 26.59
## 1  6.235                                             1 -699.943 1402.4 26.64
## 20 5.564 0.009598 0.11820                  0.003244  4 -693.567 1403.1 27.38
## 8  5.359 0.023130 0.14840 0.027410                   4 -694.253 1404.5 28.76
## 18 6.073 0.001557                          0.003082  3 -697.595 1405.2 29.44
## 5  6.207                  0.006960                   2 -699.806 1405.3 29.58
## 2  6.166 0.006571                                    2 -699.893 1405.5 29.75
## 22 5.845 0.004883         0.029770         0.004657  4 -695.778 1407.6 31.81
## 6  6.123 0.007845         0.007488                   3 -699.737 1409.5 33.72
##    weight
## 31  0.981
## 29  0.011
## 15  0.003
## 32  0.001
## 23  0.001
## 13  0.001
## 11  0.000
## 9   0.000
## 25  0.000
## 30  0.000
## 27  0.000
## 16  0.000
## 10  0.000
## 14  0.000
## 12  0.000
## 19  0.000
## 24  0.000
## 3   0.000
## 26  0.000
## 7   0.000
## 17  0.000
## 21  0.000
## 4   0.000
## 28  0.000
## 1   0.000
## 20  0.000
## 8   0.000
## 18  0.000
## 5   0.000
## 2   0.000
## 22  0.000
## 6   0.000
## Models ranked by AICc(x)
#summary(model.avg(model_dredge, subset = delta <= 2))
car::vif(get.models(model_dredge, 1)[[1]]) # looking for less than 2
##      log_tdh      log_tdr log_TotalSub     tree100m 
##     1.235699     1.878708     1.218450     1.489221
# Urban Center
global <- glm(round(RetreatDist, 0) ~ log_TotalSub + log_tdr + log_tdh + spec_rad + log_patch + road_length_m, data = near_center, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_wd_c <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_wd_c)
out
## 
## Call:
## glm(formula = round(RetreatDist, 0) ~ log_TotalSub + road_length_m + 
##     spec_rad + 1, family = poisson, data = near_center)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -4.654  -3.928  -1.932   1.507   9.858  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   25.656433   1.978299  12.969   <2e-16 ***
## log_TotalSub  -0.949802   0.081234 -11.692   <2e-16 ***
## road_length_m -0.028700   0.001686 -17.019   <2e-16 ***
## spec_rad      -0.116126   0.013624  -8.523   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1550.31  on 10  degrees of freedom
## Residual deviance:  250.73  on  7  degrees of freedom
## AIC: 323.93
## 
## Number of Fisher Scoring iterations: 5
info <- out[[12]]
#TotalSub
exp(info[2, 1])
## [1] 0.3868175
exp(info[2, 2])
## [1] 1.084625
#spec_rad
info[4, 1]
## [1] -0.1161259
info[4, 2]
## [1] 0.01362443
#road_length_m
info[3, 1]
## [1] -0.02870023
info[3, 2]
## [1] 0.001686396
with(summary(top_model_wd_c), 1 - deviance/null.deviance)
## [1] 0.8382736
model.sel(model_dredge) #estimates same signs
## Global model call: glm(formula = round(RetreatDist, 0) ~ log_TotalSub + log_tdr + 
##     log_tdh + spec_rad + log_patch + road_length_m, family = poisson, 
##     data = near_center)
## ---
## Model selection table 
##      (Int)   log_ptc   log_tdh   log_tdr log_TtS rod_lng_m spc_rad df   logLik
## 57 25.6600                               -0.9498  -0.02870 -0.1161  4 -157.967
## 58 30.1500 -0.152100                     -1.1060  -0.02902 -0.1381  5 -154.713
## 61 24.5500                      0.043540 -1.0180  -0.02943 -0.1096  5 -157.053
## 59 25.2800            0.039410           -0.9449  -0.02824 -0.1150  5 -157.751
## 62 29.2100 -0.141300            0.023170 -1.1300  -0.02941 -0.1329  6 -154.471
## 60 29.8100 -0.150300  0.032670           -1.1020  -0.02861 -0.1371  6 -154.557
## 63 23.4200            0.080480  0.059620 -1.0310  -0.02872 -0.1052  6 -156.288
## 64 28.1400 -0.132300  0.057050  0.035760 -1.1350  -0.02889 -0.1284  7 -154.072
## 31  7.1510            0.197100  0.181200 -1.3260  -0.03741          5 -182.505
## 29  8.3050                      0.158200 -1.3290  -0.04028          4 -187.159
## 30  7.1210  0.138900            0.153100 -1.1650  -0.03796          5 -183.907
## 32  5.9130  0.138500  0.205200  0.178400 -1.1510  -0.03505          6 -179.113
## 26  7.4630  0.165700                     -0.9463  -0.03684          4 -198.754
## 28  6.8810  0.162600  0.113200           -0.9302  -0.03524          5 -197.129
## 25  8.8960                               -1.1450  -0.03945          3 -203.817
## 27  8.2970            0.111800           -1.1350  -0.03776          4 -202.156
## 53 31.3700                     -0.137400          -0.02833 -0.1547  4 -212.808
## 54 28.8100  0.066910           -0.122000          -0.02833 -0.1413  5 -211.962
## 55 31.7100           -0.025670 -0.142100          -0.02855 -0.1559  5 -212.709
## 50 23.6000  0.156500                              -0.03006 -0.1142  4 -220.854
## 56 28.8400  0.066540 -0.001317 -0.122300          -0.02834 -0.1414  6 -211.961
## 52 22.1500  0.166900  0.117200                    -0.02860 -0.1091  5 -218.580
## 49 28.7400                                        -0.03090 -0.1409  3 -226.061
## 51 28.1000            0.087040                    -0.02970 -0.1397  4 -224.683
## 20  4.0230  0.396700  0.200700                    -0.03318          4 -247.643
## 24  3.6080  0.403800  0.233600  0.040720          -0.03290          5 -246.520
## 18  5.0340  0.402000                              -0.03605          3 -253.698
## 22  5.0100  0.402600            0.004235          -0.03608          4 -253.682
## 44 59.3200 -0.294400  0.290300           -1.2290           -0.3456  5 -280.443
## 48 57.9200 -0.274000  0.318300  0.034050 -1.2420           -0.3391  6 -280.067
## 19  7.6510            0.117000                    -0.04320          3 -296.354
## 17  8.2780                                        -0.04510          2 -298.758
## 47 45.6500            0.404600  0.120100 -1.0110           -0.2768  5 -291.713
## 42 65.3400 -0.286500                     -1.2440           -0.3775  4 -295.598
## 23  7.4970            0.128500  0.021760          -0.04330          4 -295.974
## 21  8.2470                      0.007374          -0.04520          3 -298.710
## 46 67.0600 -0.327400           -0.064060 -1.2100           -0.3850  5 -293.606
## 43 47.6400            0.302300           -0.8809           -0.2838  4 -297.299
## 41 54.4600                               -0.9217           -0.3201  3 -311.598
## 45 54.4700                      0.011800 -0.9369           -0.3204  4 -311.519
## 39 51.1800            0.254800 -0.088140                   -0.3079  4 -369.806
## 35 50.2900            0.332300                             -0.3072  3 -374.060
## 36 47.2000  0.061770  0.340300                             -0.2898  4 -372.810
## 40 51.5300 -0.006439  0.251600 -0.091010                   -0.3097  5 -369.797
## 37 56.2400                     -0.148700                   -0.3321  3 -379.890
## 38 59.4600 -0.071570           -0.171900                   -0.3494  4 -378.615
## 33 57.0600                                                 -0.3426  2 -397.328
## 34 54.3200  0.060710                                       -0.3272  3 -396.159
## 16 -3.9730  0.508400  1.009000  0.383500 -0.7939                    5 -460.224
## 12 -0.3822  0.396200  0.739500           -0.5345                    4 -519.353
## 8  -4.2820  0.633100  0.895500  0.241500                            4 -519.762
## 15  0.8805            0.873600  0.262600 -1.2670                    4 -546.164
## 4  -1.5460  0.523300  0.707700                                      3 -549.326
## 11  2.6130            0.727600           -1.0120                    3 -581.675
## 14  1.3650  0.495600            0.134500 -0.5765                    4 -605.004
## 10  2.1910  0.458700                     -0.4552                    3 -617.478
## 2   1.0670  0.569500                                                2 -644.911
## 6   0.7389  0.588600            0.039460                            3 -643.523
## 13  5.1230                      0.088690 -0.9794                    3 -686.072
## 9   5.5000                               -0.8882                    2 -692.365
## 3   2.3230            0.573700                                      2 -717.922
## 7   2.2140            0.580600  0.015630                            3 -717.721
## 5   5.0070                     -0.063300                            2 -803.362
## 1   4.6770                                                          1 -807.758
##      AICc   delta weight
## 57  330.6    0.00  0.567
## 58  331.4    0.83  0.375
## 61  336.1    5.51  0.036
## 59  337.5    6.90  0.018
## 62  341.9   11.34  0.002
## 60  342.1   11.51  0.002
## 63  345.6   14.98  0.000
## 64  359.5   28.88  0.000
## 31  387.0   56.41  0.000
## 29  389.0   58.38  0.000
## 30  389.8   59.21  0.000
## 32  391.2   60.63  0.000
## 26  412.2   81.58  0.000
## 28  416.3   85.66  0.000
## 25  417.1   86.46  0.000
## 27  419.0   88.38  0.000
## 53  440.3  109.68  0.000
## 54  445.9  115.32  0.000
## 55  447.4  116.82  0.000
## 50  456.4  125.77  0.000
## 56  456.9  126.32  0.000
## 52  459.2  128.56  0.000
## 49  461.6  130.95  0.000
## 51  464.0  133.43  0.000
## 20  510.0  179.35  0.000
## 24  515.0  184.44  0.000
## 18  516.8  186.22  0.000
## 22  522.0  191.43  0.000
## 44  582.9  252.28  0.000
## 48  593.1  262.53  0.000
## 19  602.1  271.54  0.000
## 17  603.0  272.42  0.000
## 47  605.4  274.82  0.000
## 42  605.9  275.26  0.000
## 23  606.6  276.01  0.000
## 21  606.8  276.25  0.000
## 46  609.2  278.61  0.000
## 43  609.3  278.66  0.000
## 41  632.6  302.02  0.000
## 45  637.7  307.10  0.000
## 39  754.3  423.68  0.000
## 35  757.5  426.95  0.000
## 36  760.3  429.69  0.000
## 40  761.6  430.99  0.000
## 37  769.2  438.61  0.000
## 38  771.9  441.30  0.000
## 33  800.2  469.56  0.000
## 34  801.7  471.15  0.000
## 16  942.4  611.85  0.000
## 12 1053.4  722.77  0.000
## 8  1054.2  723.59  0.000
## 15 1107.0  776.39  0.000
## 4  1108.1  777.48  0.000
## 11 1172.8  842.18  0.000
## 14 1224.7  894.07  0.000
## 10 1244.4  913.78  0.000
## 2  1295.3  964.72  0.000
## 6  1296.5  965.87  0.000
## 13 1381.6 1050.97  0.000
## 9  1390.2 1059.63  0.000
## 3  1441.3 1110.74  0.000
## 7  1444.9 1114.27  0.000
## 5  1612.2 1281.62  0.000
## 1  1618.0 1287.36  0.000
## Models ranked by AICc(x)
out <- summary(model.avg(model_dredge, subset = delta <= 2))
out
## 
## Call:
## model.avg(object = model_dredge, subset = delta <= 2)
## 
## Component model call: 
## glm(formula = round(RetreatDist, 0) ~ <2 unique rhs>, family = poisson, 
##      data = near_center)
## 
## Component models: 
##      df  logLik   AICc delta weight
## 234   4 -157.97 330.60  0.00    0.6
## 1234  5 -154.71 331.43  0.83    0.4
## 
## Term codes: 
##     log_patch  log_TotalSub road_length_m      spec_rad 
##             1             2             3             4 
## 
## Model-averaged coefficients:  
## (full average) 
##                Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   27.444071   3.180180    3.580799   7.664   <2e-16 ***
## log_TotalSub  -1.011887   0.119466    0.136311   7.423   <2e-16 ***
## road_length_m -0.028829   0.001723    0.002106  13.691   <2e-16 ***
## spec_rad      -0.124872   0.018357    0.021193   5.892   <2e-16 ***
## log_patch     -0.060550   0.083552    0.088232   0.686    0.493    
##  
## (conditional average) 
##                Estimate Std. Error Adjusted SE z value Pr(>|z|)    
## (Intercept)   27.444071   3.180180    3.580799   7.664   <2e-16 ***
## log_TotalSub  -1.011887   0.119466    0.136311   7.423   <2e-16 ***
## road_length_m -0.028829   0.001723    0.002106  13.691   <2e-16 ***
## spec_rad      -0.124872   0.018357    0.021193   5.892   <2e-16 ***
## log_patch     -0.152072   0.060123    0.075060   2.026   0.0428 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
info <- out[[10]]
#TotalSub
exp(info[2, 1])
## [1] 0.3635323
exp(info[2, 2])
## [1] 1.126895
#spec_rad
info[4, 1]
## [1] -0.1248715
info[4, 2]
## [1] 0.0183575
#road_length_m
info[3, 1]
## [1] -0.02882948
info[3, 2]
## [1] 0.001722559
#log_patch
exp(info[5, 1])
## [1] 0.8589265
exp(info[5, 2])
## [1] 1.061967
car::vif(get.models(model_dredge, 1)[[1]]) # looking for less than 2
##  log_TotalSub road_length_m      spec_rad 
##      1.046036      1.945607      2.004812
# If any value is much greater than 2, this means there is quite a bit of multicollinearity, and we can try removing the variable with the highest vif. Here is the code you would run to accomplish this.

#global2 <- update(global, .~. -spec_rad)
#options(na.action = "na.fail") # Required for dredge to run
#model_dredge <- dredge(global2, beta = "none", evaluate = T, rank = AICc)
#options(na.action = "na.omit") # set back to default
#top_model_wd_c <- get.models(model_dredge, subset = 1)[[1]]
#out <- summary(top_model_wd_c)
#out
#info <- out[[12]]
#car::vif(get.models(model_dredge, 1)[[1]]) 

Overall Coefficients: Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.38311 0.06200 86.83 <2e-16 log_tdr -0.18279 0.01024 -17.84 <2e-16 negative log_TotalSub 0.57484 0.01653 34.78 <2e-16 *** positive

R-squared = 0.3424213

Urban Forest Coefficients: Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.247233 0.253935 20.664 < 2e-16 log_tdh 0.192623 0.045552 4.229 2.35e-05 positive log_tdr 0.097475 0.018717 5.208 1.91e-07 *** positive log_TotalSub -0.182147 0.037627 -4.841 1.29e-06 *** negative tree100m 0.007928 0.001747 4.540 5.64e-06 *** positive

R-squared = 0.03720731

Urban Center Coefficients: Estimate Std. Error z value Pr(>|z|)
(Intercept) 25.656433 1.978299 12.969 <2e-16 log_TotalSub -0.949802 0.081234 -11.692 <2e-16 negative road_length_m -0.028700 0.001686 -17.019 <2e-16 *** negative spec_rad -0.116126 0.013624 -8.523 <2e-16 *** negative

R-squared = 0.8382736

(conditional average) Estimate Std. Error Adjusted SE z value Pr(>|z|)
(Intercept) 27.444071 3.180180 3.580799 7.664 <2e-16 log_TotalSub -1.011887 0.119466 0.136311 7.423 <2e-16 road_length_m -0.028829 0.001723 0.002106 13.691 <2e-16 spec_rad -0.124872 0.018357 0.021193 5.892 <2e-16 log_patch -0.152072 0.060123 0.075060 2.026 0.0428 *

Web Distance by Land

dist_land <- glm(round(RetreatDist, 0) ~ Land, data = near, family = poisson)
summary(dist_land)
## 
## Call:
## glm(formula = round(RetreatDist, 0) ~ Land, family = poisson, 
##     data = near)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -15.5268   -8.9632   -0.6712    3.7733   24.6953  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           4.12713    0.06350  64.994  < 2e-16 ***
## LandUrban, Medium     1.05689    0.07179  14.723  < 2e-16 ***
## LandUrban, Low       -1.08261    0.16686  -6.488 8.69e-11 ***
## LandDeciduous Forest  2.19774    0.06548  33.561  < 2e-16 ***
## LandWoody Wetlands    1.86015    0.06978  26.658  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 5895.3  on 20  degrees of freedom
## Residual deviance: 2266.0  on 16  degrees of freedom
## AIC: 2420.5
## 
## Number of Fisher Scoring iterations: 5
car::Anova(dist_land, test.statistic = "LR")
## Analysis of Deviance Table (Type II tests)
## 
## Response: round(RetreatDist, 0)
##      LR Chisq Df Pr(>Chisq)    
## Land   3629.3  4  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
comp <- glht(dist_land, linfct = mcp(Land = "Tukey"))
summary(comp)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: glm(formula = round(RetreatDist, 0) ~ Land, family = poisson, 
##     data = near)
## 
## Linear Hypotheses:
##                                        Estimate Std. Error z value Pr(>|z|)    
## Urban, Medium - Urban, High == 0        1.05689    0.07179  14.723  < 1e-10 ***
## Urban, Low - Urban, High == 0          -1.08261    0.16686  -6.488  3.9e-10 ***
## Deciduous Forest - Urban, High == 0     2.19774    0.06548  33.561  < 1e-10 ***
## Woody Wetlands - Urban, High == 0       1.86015    0.06978  26.658  < 1e-10 ***
## Urban, Low - Urban, Medium == 0        -2.13951    0.15789 -13.550  < 1e-10 ***
## Deciduous Forest - Urban, Medium == 0   1.14084    0.03711  30.745  < 1e-10 ***
## Woody Wetlands - Urban, Medium == 0     0.80326    0.04425  18.154  < 1e-10 ***
## Deciduous Forest - Urban, Low == 0      3.28035    0.15513  21.146  < 1e-10 ***
## Woody Wetlands - Urban, Low == 0        2.94277    0.15699  18.745  < 1e-10 ***
## Woody Wetlands - Deciduous Forest == 0 -0.33758    0.03306 -10.212  < 1e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
near <- near %>% 
  mutate(Land = fct_relevel(Land, "Deciduous Forest", 
                            "Woody Wetlands", 
                            "Urban, High", 
                            "Urban, Medium", 
                            "Urban, Low"))
nd <- data.frame(Land = factor(levels(sites$Land)))
predictions <- augment(dist_land, newdata = nd, se_fit = TRUE, type.predict = "response")
predictions <- predictions %>% 
  rename("response" = ".fitted", "SE" = ".se.fit")
#predictions <- summary(emmeans(dist_land, ~Land),type = "response")

ggplot(near, aes(x = Land, y = RetreatDist)) + 
  geom_jitter(color = "grey", width = 0.1, size = 0.5) +
  geom_point(aes(x = Land, y = response, color = Land), size = 1, data = predictions) + 
  geom_errorbar(aes(x = Land, 
                  ymin = response - SE, 
                  ymax = response + SE, color = Land), data = predictions, inherit.aes = FALSE, width = 0.25, size = 1) +
  theme_classic() +
  scale_color_manual(values = c("Deciduous Forest" = "#1b9e77", "Woody Wetlands" = "#1b9e77", "Urban, High" = "#d95f02", "Urban, Medium" = "#d95f02", "Urban, Low" = "#d95f02")) +
  xlab("Land Cover Class") +
  ylab("Distance from Focal Web \nto Neighbor Web (cm)") + 
  theme(text = element_text(size = 10)) + 
  theme(axis.text.x=element_text(colour="black", size=10)) + 
  theme(axis.text.y=element_text(colour="black", size=10)) +
  theme(legend.position = "none") +
  theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust=1)) +
  annotate(geom="text", x=1, y=max(sites$NumSpiders), label="A", color="black", size = 3) +
  annotate(geom="text", x=2, y=max(sites$NumSpiders), label="B", color="black", size = 3) +
  annotate(geom="text", x=3, y=max(sites$NumSpiders), label="C", color="black", size = 3) +
  annotate(geom="text", x=4, y=max(sites$NumSpiders), label="D", color="black", size = 3) +
  annotate(geom="text", x=5, y=max(sites$NumSpiders), label="E", color="black", size = 3)

Distance between webs significantly varies by land cover class (Chi = 3629.3, df = 4, p < 0.001).

Response Variable: Web Height

Web Height by Location

webs$RetreatHeight[webs$RetreatHeight == 0] <- 0.1
ggarrange(ggplot(webs, aes(x = RetreatHeight)) +
            geom_density(),
          ggplot(webs, aes(x = log(RetreatHeight))) +
            geom_density(), ncol = 2)

shapiro.test(webs$RetreatHeight) # not normal
## 
##  Shapiro-Wilk normality test
## 
## data:  webs$RetreatHeight
## W = 0.7433, p-value = 7.58e-14
shapiro.test(log(webs$RetreatHeight)) # not normal
## 
##  Shapiro-Wilk normality test
## 
## data:  log(webs$RetreatHeight)
## W = 0.96991, p-value = 0.005222
height <- glmer(round(RetreatHeight, 0) ~ Location + (1 | ID), data = webs, family = poisson)
summary(height)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: round(RetreatHeight, 0) ~ Location + (1 | ID)
##    Data: webs
## 
##      AIC      BIC   logLik deviance df.resid 
##   2455.2   2463.8  -1224.6   2449.2      128 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.4663 -2.2974 -0.5413  1.2882 23.4563 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 0.5289   0.7273  
## Number of obs: 131, groups:  ID, 22
## 
## Fixed effects:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            3.5235     0.2331  15.115   <2e-16 ***
## LocationUrban Center  -0.7527     0.3157  -2.384   0.0171 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## LctnUrbnCnt -0.738
height2 <- lmer(log(RetreatHeight) ~ Location + (1 | ID), data = webs)
summary(height2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: log(RetreatHeight) ~ Location + (1 | ID)
##    Data: webs
## 
## REML criterion at convergence: 369.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2771 -0.5190  0.0362  0.5077  3.1262 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.4884   0.6988  
##  Residual             0.7789   0.8826  
## Number of obs: 131, groups:  ID, 22
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)            3.3407     0.2750  12.150
## LocationUrban Center  -0.8710     0.3567  -2.442
## 
## Correlation of Fixed Effects:
##             (Intr)
## LctnUrbnCnt -0.771

Webs are significantly higher in the urban forest than the urban center (z = -2.384, df = 1, 21, p = 0.017).

Web Height by Predictors

#overall
models = list(ceiling(RetreatHeight) ~ 1,
              ceiling(RetreatHeight) ~ log_TotalSub,
              ceiling(RetreatHeight) ~ log_tdr, 
              ceiling(RetreatHeight) ~ log_tdh,
              ceiling(RetreatHeight) ~ log_TotalSub + log_tdr,
              ceiling(RetreatHeight) ~ log_TotalSub + log_tdh,
              ceiling(RetreatHeight) ~ log_tdr + log_tdh,
              ceiling(RetreatHeight) ~ log_TotalSub + log_tdr + log_tdh)
fits = lapply(models, glm, data = webs, family = "poisson")
modnames = sapply(models, function(ff)deparse(ff[[3]]))
pander(aictab(fits, modname = modnames), caption="Table 1. AICc model selection table for Search Distance to the Focal Web.", split.tables = Inf)
Table 1. AICc model selection table for Search Distance to the Focal Web.
  Modnames K AICc Delta_AICc ModelLik AICcWt LL Cum.Wt
8 log_TotalSub + log_tdr + log_tdh 4 3729 0 1 0.9998 -1860 0.9998
5 log_TotalSub + log_tdr 3 3746 17.36 0.0001696 0.0001696 -1870 1
6 log_TotalSub + log_tdh 3 3864 135.5 3.853e-30 3.852e-30 -1929 1
2 log_TotalSub 2 3950 220.8 1.139e-48 1.139e-48 -1973 1
7 log_tdr + log_tdh 3 4098 369 7.524e-81 7.523e-81 -2046 1
3 log_tdr 2 4114 385 2.519e-84 2.518e-84 -2055 1
1 1 1 4436 707.4 2.468e-154 2.467e-154 -2217 1
4 log_tdh 2 4438 709.1 1.039e-154 1.039e-154 -2217 1
global <- glmer(round(RetreatHeight, 0) ~ log_TotalSub + log_tdr + log_tdh + (1 | ID), data = webs, family = poisson)
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = TRUE, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_wh_o <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_wh_o)
out
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: round(RetreatHeight, 0) ~ log_tdr + log_TotalSub + (1 | ID)
##    Data: webs
## 
##      AIC      BIC   logLik deviance df.resid 
##   2451.0   2462.5  -1221.5   2443.0      127 
## 
## Scaled residuals: 
##    Min     1Q Median     3Q    Max 
## -7.455 -2.323 -0.569  1.283 23.499 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 0.3959   0.6292  
## Number of obs: 131, groups:  ID, 22
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    3.1686     0.6215   5.098 3.43e-07 ***
## log_tdr       -0.1896     0.1023  -1.854  0.06378 .  
## log_TotalSub   0.4812     0.1659   2.901  0.00372 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_tdr
## log_tdr     -0.864       
## log_TotalSb -0.621  0.205
info <- out[[10]]
# TotalSub
exp(info[3, 1])
## [1] 1.61805
exp(info[3, 2])
## [1] 1.180415
#Traffic_Dist_Road
exp(info[2, 1])
## [1] 0.8272671
exp(info[2, 1])
## [1] 0.8272671
r.squaredGLMM(top_model_wh_o) 
##                 R2m       R2c
## delta     0.3674559 0.9528885
## lognormal 0.3677262 0.9535893
## trigamma  0.3671774 0.9521663
model.sel(model_dredge) #estimates same signs
## Global model call: glmer(formula = round(RetreatHeight, 0) ~ log_TotalSub + log_tdr + 
##     log_tdh + (1 | ID), data = webs, family = poisson)
## ---
## Model selection table 
##   (Int) log_tdh log_tdr log_TtS df    logLik   AICc delta weight
## 7 3.169         -0.1896  0.4812  4 -1221.517 2451.4  0.00  0.340
## 8 4.672 -0.3307 -0.2116  0.3900  5 -1220.799 2452.1  0.73  0.236
## 5 2.171                  0.5444  3 -1223.114 2452.4  1.07  0.199
## 6 3.135 -0.2298          0.4860  4 -1222.810 2453.9  2.59  0.093
## 4 6.566 -0.5842 -0.2691          4 -1222.978 2454.3  2.92  0.079
## 3 4.285         -0.2509          3 -1225.073 2456.3  4.98  0.028
## 2 5.108 -0.5315                  3 -1225.720 2457.6  6.28  0.015
## 1 3.110                          2 -1227.125 2458.3  6.99  0.010
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | ID
out <- summary(model.avg(model_dredge, subset = delta <= 2))
out
## 
## Call:
## model.avg(object = model_dredge, subset = delta <= 2)
## 
## Component model call: 
## glmer(formula = round(RetreatHeight, 0) ~ <3 unique rhs>, data = webs, 
##      family = poisson)
## 
## Component models: 
##     df   logLik    AICc delta weight
## 23   4 -1221.52 2451.35  0.00   0.44
## 123  5 -1220.80 2452.08  0.73   0.30
## 3    3 -1223.11 2452.42  1.07   0.26
## 
## Term codes: 
##      log_tdh      log_tdr log_TotalSub 
##            1            2            3 
## 
## Model-averaged coefficients:  
## (full average) 
##              Estimate Std. Error Adjusted SE z value Pr(>|z|)   
## (Intercept)    3.3701     1.2954      1.3012   2.590  0.00959 **
## log_tdr       -0.1476     0.1237      0.1243   1.187  0.23518   
## log_TotalSub   0.4697     0.1814      0.1830   2.567  0.01026 * 
## log_tdh       -0.1008     0.2138      0.2148   0.469  0.63894   
##  
## (conditional average) 
##              Estimate Std. Error Adjusted SE z value Pr(>|z|)   
## (Intercept)    3.3701     1.2954      1.3012   2.590  0.00959 **
## log_tdr       -0.1986     0.1022      0.1032   1.924  0.05432 . 
## log_TotalSub   0.4697     0.1814      0.1830   2.567  0.01026 * 
## log_tdh       -0.3307     0.2719      0.2745   1.205  0.22832   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
info <- out[[10]]
#tdr
exp(info[2, 1])
## [1] 0.8198516
exp(info[2, 2])
## [1] 1.107659
#TotalSub
exp(info[3, 1])
## [1] 1.599464
exp(info[3, 2])
## [1] 1.198902
#tdh
exp(info[4, 1])
## [1] 0.7184019
exp(info[4, 2])
## [1] 1.312448
car::vif(get.models(model_dredge, 5)[[1]]) # looking for less than 2
##  log_tdh  log_tdr 
## 1.006424 1.006424
# Urban Forest
webs_forest <- webs_forest %>% 
  mutate(tree100m_frac = I(tree100m/100))
global <- glmer(round(RetreatHeight, 0) ~ tree100m_frac + log_TotalSub + log_tdr + log_tdh + log_patch + (1 | ID), data = webs_forest, family = poisson, control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e4)))
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = TRUE, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_wh_f <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_wh_f)
out
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: round(RetreatHeight, 0) ~ tree100m_frac + (1 | ID)
##    Data: webs_forest
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000))
## 
##      AIC      BIC   logLik deviance df.resid 
##    776.9    781.1   -385.5    770.9       27 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.4737 -2.5373 -0.4187  1.6654 14.8549 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 0.2732   0.5227  
## Number of obs: 30, groups:  ID, 10
## 
## Fixed effects:
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept)     1.1283     0.8355   1.350  0.17687   
## tree100m_frac   5.0864     1.7266   2.946  0.00322 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## tre100m_frc -0.979
info <- out[[10]]

#tree
exp(info[2, 1])
## [1] 161.8125
exp(info[2, 2])
## [1] 5.621704
r.squaredGLMM(top_model_wh_f) 
##                 R2m       R2c
## delta     0.4789921 0.9601273
## lognormal 0.4792051 0.9605545
## trigamma  0.4787743 0.9596909
model.sel(model_dredge) #estimates same signs
## Global model call: glmer(formula = round(RetreatHeight, 0) ~ tree100m_frac + log_TotalSub + 
##     log_tdr + log_tdh + log_patch + (1 | ID), data = webs_forest, 
##     family = poisson, control = glmerControl(optimizer = "bobyqa", 
##         optCtrl = list(maxfun = 20000)))
## ---
## Model selection table 
##      (Int)   log_ptc   log_tdh  log_tdr log_TtS t10_frc df   logLik  AICc delta
## 17  1.1280                                        5.086  3 -385.466 777.9  0.00
## 25 -0.5283                               0.6026   5.398  4 -384.247 778.1  0.24
## 29  0.4737                     -0.24150  0.8134   4.187  5 -383.312 779.1  1.27
## 21  1.6770                     -0.08359           4.627  4 -385.368 780.3  2.48
## 18  0.6326  0.048830                              5.048  4 -385.447 780.5  2.64
## 19  1.4150           -0.081090                    5.082  4 -385.453 780.5  2.65
## 26 -0.9924  0.045900                     0.6020   5.361  5 -384.225 781.0  3.09
## 27 -0.7050            0.045410           0.6084   5.403  5 -384.241 781.0  3.13
## 13  3.0190                     -0.46340  0.9361          4 -385.833 781.3  3.41
## 1   3.5230                                               2 -388.596 781.6  3.78
## 31  1.8200           -0.286300 -0.29860  0.8278   3.868  6 -383.114 781.9  4.02
## 5   4.7410                     -0.30630                  3 -387.544 782.0  4.16
## 30  0.2346  0.022980           -0.23990  0.8118   4.176  6 -383.306 782.3  4.41
## 15  5.4390           -0.598600 -0.54560  0.9406          5 -385.233 783.0  5.11
## 23  2.8310           -0.241900 -0.12900           4.366  5 -385.272 783.0  5.19
## 22  1.2430  0.041400           -0.08147           4.606  5 -385.354 783.2  5.35
## 20  0.9268  0.043700 -0.068410                    5.048  5 -385.438 783.4  5.52
## 9   2.3670                               0.4615          3 -388.243 783.4  5.55
## 7   7.1380           -0.591800 -0.38600                  4 -387.126 783.9  6.00
## 2   2.1190  0.133400                                     3 -388.521 784.0  6.11
## 3   3.9700           -0.127100                           3 -388.578 784.1  6.22
## 28 -1.2720  0.050350  0.060320           0.6096   5.364  6 -384.216 784.1  6.23
## 14  2.4630  0.052260           -0.45870  0.9313          5 -385.814 784.1  6.27
## 6   3.8740  0.080080           -0.30040                  4 -387.511 784.6  6.77
## 32  1.8690 -0.003859 -0.288200 -0.29930  0.8281   3.868  7 -383.114 785.3  7.46
## 10  0.9673  0.133100                     0.4611          4 -388.163 785.9  8.07
## 11  2.4950           -0.033120           0.4569          4 -388.242 786.1  8.23
## 16  5.5470 -0.008644 -0.602600 -0.54690  0.9414          6 -385.233 786.1  8.26
## 24  2.5780  0.019930 -0.232700 -0.12620           4.366  6 -385.269 786.2  8.33
## 4   2.5080  0.126600 -0.089910                           4 -388.512 786.6  8.77
## 8   6.8700  0.021180 -0.582100 -0.38310                  5 -387.123 786.7  8.89
## 12  0.9340  0.133700  0.007122           0.4621          5 -388.163 788.8 10.97
##    weight
## 17  0.196
## 25  0.174
## 29  0.104
## 21  0.057
## 18  0.052
## 19  0.052
## 26  0.042
## 27  0.041
## 13  0.036
## 1   0.030
## 31  0.026
## 5   0.025
## 30  0.022
## 15  0.015
## 23  0.015
## 22  0.014
## 20  0.012
## 9   0.012
## 7   0.010
## 2   0.009
## 3   0.009
## 28  0.009
## 14  0.009
## 6   0.007
## 32  0.005
## 10  0.003
## 11  0.003
## 16  0.003
## 24  0.003
## 4   0.002
## 8   0.002
## 12  0.001
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | ID
out <- summary(model.avg(model_dredge, subset = delta <= 2))
out
## 
## Call:
## model.avg(object = model_dredge, subset = delta <= 2)
## 
## Component model call: 
## glmer(formula = round(RetreatHeight, 0) ~ <3 unique rhs>, data = 
##      webs_forest, family = poisson, control = glmerControl(optimizer = 
##      "bobyqa", optCtrl = list(maxfun = 20000)))
## 
## Component models: 
##     df  logLik   AICc delta weight
## 3    3 -385.47 777.86  0.00   0.41
## 23   4 -384.25 778.09  0.24   0.37
## 123  5 -383.31 779.12  1.27   0.22
## 
## Term codes: 
##       log_tdr  log_TotalSub tree100m_frac 
##             1             2             3 
## 
## Model-averaged coefficients:  
## (full average) 
##               Estimate Std. Error Adjusted SE z value Pr(>|z|)   
## (Intercept)    0.37641    1.34169     1.38813   0.271    0.786   
## tree100m_frac  5.00334    1.70559     1.78233   2.807    0.005 **
## log_TotalSub   0.39968    0.44424     0.45314   0.882    0.378   
## log_tdr       -0.05297    0.12729     0.12981   0.408    0.683   
##  
## (conditional average) 
##               Estimate Std. Error Adjusted SE z value Pr(>|z|)   
## (Intercept)     0.3764     1.3417      1.3881   0.271   0.7863   
## tree100m_frac   5.0033     1.7056      1.7823   2.807   0.0050 **
## log_TotalSub    0.6815     0.3801      0.3976   1.714   0.0865 . 
## log_tdr        -0.2415     0.1683      0.1769   1.365   0.1721   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
info <- out[[10]]
#tree
exp(info[2, 1])
## [1] 148.9091
exp(info[2, 2])
## [1] 5.504634
#TotalSub
exp(info[3, 1])
## [1] 1.976754
exp(info[3, 2])
## [1] 1.462413
#tdr
exp(info[4, 1])
## [1] 0.7854307
exp(info[4, 2])
## [1] 1.183328
#car::vif(get.models(model_dredge, 1)[[1]]) # looking for less than 2

# Urban Center
global <- glmer(round(RetreatHeight, 0) ~ log_TotalSub + log_tdr + log_tdh + scale(spec_rad, center = FALSE) + log_patch + scale(road_length_m, center = FALSE) + (1 | ID), data = webs_center, family = poisson,  control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e4)))
options(na.action = "na.fail") # Required for dredge to run
model_dredge <- dredge(global, beta = "none", evaluate = T, rank = AICc)
options(na.action = "na.omit") # set back to default
top_model_wh_c <- get.models(model_dredge, subset = 1)[[1]]
out <- summary(top_model_wh_c)
out
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: round(RetreatHeight, 0) ~ log_TotalSub + (1 | ID)
##    Data: webs_center
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1675.6   1683.4   -834.8   1669.6       98 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.9527 -2.1754 -0.7688  1.1157 23.3877 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 0.427    0.6534  
## Number of obs: 101, groups:  ID, 12
## 
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    2.0705     0.4687   4.417    1e-05 ***
## log_TotalSub   0.6449     0.3921   1.645      0.1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## log_TotalSb -0.912
info <- out[[10]]
#TotalSub
exp(info[2, 1])
## [1] 1.905754
exp(info[2, 2])
## [1] 1.480087
r.squaredGLMM(top_model_wh_c) 
##                 R2m       R2c
## delta     0.1335054 0.9123246
## lognormal 0.1337785 0.9141914
## trigamma  0.1332201 0.9103755
model.sel(model_dredge) #estimates same signs
## Global model call: glmer(formula = round(RetreatHeight, 0) ~ log_TotalSub + log_tdr + 
##     log_tdh + scale(spec_rad, center = FALSE) + log_patch + scale(road_length_m, 
##     center = FALSE) + (1 | ID), data = webs_center, family = poisson, 
##     control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 20000)))
## ---
## Model selection table 
##      (Int)   log_ptc log_tdh  log_tdr log_TtS scl(rod_lng_m,F) scl(spc_rad,F)
## 9   2.0700                             0.6449                                
## 17  1.5220                                              1.2640               
## 25  1.1270                             0.5495           1.0610               
## 1   2.7710                                                                   
## 15  5.0850           -0.4666 -0.27140  0.8864                                
## 13  2.9020                   -0.20300  0.8655                                
## 29  1.9570                   -0.18140  0.7567           0.9623               
## 10  0.5941  0.212700                   0.8131                                
## 11  3.2880           -0.2983           0.6102                                
## 31  3.9960           -0.3908 -0.24340  0.7994           0.7472               
## 26 -0.2737  0.204500                   0.7136           1.0400               
## 3   4.1520           -0.3491                                                 
## 16  3.6640  0.189400 -0.4466 -0.26490  1.0320                                
## 14  1.4490  0.207200         -0.19940  1.0260                                
## 19  2.6380           -0.2503                            1.1370               
## 47 10.9400           -0.4888 -0.30510  0.8550                         -5.6100
## 27  2.1260           -0.2220           0.5334           0.9563               
## 30  0.5680  0.200600         -0.17820  0.9144           0.9438               
## 41  4.0170                             0.6217                         -1.9400
## 21  1.8900                   -0.06810                   1.2550               
## 33  7.1820                                                            -4.4560
## 5   3.1630                   -0.07440                                        
## 45  7.6020                   -0.22770  0.8381                         -4.5850
## 18  1.1550  0.058600                                    1.2750               
## 49  2.7180                                              1.2260        -1.1700
## 2   2.5110  0.042650                                                         
## 57  0.4364                             0.5557           1.0810         0.6708
## 32  2.5950  0.188000 -0.3724 -0.23720  0.9446           0.7430               
## 12  1.8130  0.201700 -0.2798           0.7718                                
## 7   5.2570           -0.4454 -0.13740                                        
## 61  4.1320                   -0.19400  0.7521           0.8959        -2.0580
## 58 -6.8380  0.290000                   0.8359           1.1990         5.8110
## 28  0.6993  0.197600 -0.2057           0.6929           0.9442               
## 35  8.4110           -0.3455                                          -4.3160
## 42 -1.8130  0.245300                   0.8647                          2.1740
## 23  3.6580           -0.3378 -0.11660                   1.0770               
## 63  8.1140           -0.4197 -0.27120  0.7946           0.6106        -3.7540
## 43  5.2640           -0.2986           0.5868                         -1.9680
## 4   3.9100  0.039120 -0.3480                                                 
## 48  7.0890  0.144400 -0.4632 -0.28430  0.9808                         -2.9570
## 46  2.7650  0.189700         -0.20600  1.0050                         -1.1630
## 37  9.4750                   -0.11310                                 -6.1690
## 20  2.2840  0.054810 -0.2479                            1.1480               
## 51  4.1220           -0.2533                            1.0890        -1.4390
## 62 -2.4800  0.239600         -0.16130  0.9506           1.0250         2.6290
## 53  4.6640                   -0.08502                   1.1680        -2.6250
## 59  1.7370           -0.2211           0.5369           0.9676         0.3741
## 39 12.6600           -0.4740 -0.18600                                 -7.1030
## 22  1.5900  0.042940         -0.06241                   1.2640               
## 34  7.2260 -0.002300                                                  -4.4850
## 6   2.9940  0.024930         -0.07111                                        
## 50  1.5630  0.054680                                    1.2620        -0.3754
## 44 -0.2937  0.230100 -0.2769           0.8170                          1.8920
## 64  2.3930  0.190500 -0.3708 -0.23590  0.9467           0.7490         0.1666
## 60 -5.4110  0.276600 -0.1858           0.8069           1.1000         5.3250
## 8   5.2250  0.004296 -0.4449 -0.13670                                        
## 55  8.1600           -0.3697 -0.14770                   0.9246        -4.1010
## 36  8.4990 -0.004638 -0.3456                                          -4.3750
## 24  3.4590  0.025370 -0.3338 -0.11270                   1.0850               
## 38 10.9000 -0.060800         -0.12730                                 -7.1610
## 52  3.1130  0.046910 -0.2498                            1.1220        -0.7542
## 40 15.2300 -0.103400 -0.4945 -0.21340                                 -8.8370
## 54  4.3050  0.013360         -0.08163                   1.1790        -2.3730
## 56  9.2880 -0.037650 -0.3816 -0.15930                   0.8855        -4.8610
##    df   logLik   AICc delta weight
## 9   3 -834.798 1675.8  0.00  0.053
## 17  3 -834.876 1676.0  0.16  0.049
## 25  4 -833.864 1676.1  0.30  0.046
## 1   2 -836.013 1676.1  0.30  0.046
## 15  5 -832.761 1676.2  0.31  0.046
## 13  4 -833.928 1676.3  0.43  0.043
## 29  5 -833.059 1676.7  0.91  0.034
## 10  4 -834.313 1677.0  1.20  0.029
## 11  4 -834.368 1677.2  1.31  0.028
## 31  6 -832.182 1677.3  1.41  0.026
## 26  5 -833.334 1677.3  1.46  0.026
## 3   3 -835.526 1677.3  1.46  0.026
## 16  6 -832.226 1677.3  1.50  0.025
## 14  5 -833.394 1677.4  1.58  0.024
## 19  4 -834.592 1677.6  1.76  0.022
## 47  6 -832.378 1677.6  1.80  0.022
## 27  5 -833.598 1677.8  1.98  0.020
## 30  6 -832.473 1677.8  2.00  0.020
## 41  4 -834.764 1677.9  2.10  0.019
## 21  4 -834.769 1678.0  2.11  0.019
## 33  3 -835.855 1678.0  2.11  0.018
## 5   3 -835.907 1678.1  2.22  0.018
## 45  5 -833.716 1678.1  2.22  0.018
## 18  4 -834.835 1678.1  2.24  0.017
## 49  4 -834.864 1678.1  2.30  0.017
## 2   3 -835.995 1678.2  2.39  0.016
## 57  5 -833.859 1678.3  2.51  0.015
## 32  7 -831.591 1678.4  2.54  0.015
## 12  5 -833.904 1678.4  2.59  0.015
## 7   4 -835.167 1678.8  2.91  0.012
## 61  6 -833.014 1678.9  3.08  0.011
## 58  6 -833.052 1679.0  3.15  0.011
## 28  6 -833.085 1679.1  3.22  0.011
## 35  4 -835.366 1679.1  3.30  0.010
## 42  5 -834.278 1679.2  3.34  0.010
## 23  5 -834.293 1679.2  3.37  0.010
## 63  7 -832.016 1679.2  3.39  0.010
## 43  5 -834.330 1679.3  3.45  0.009
## 4   4 -835.510 1679.4  3.59  0.009
## 48  7 -832.144 1679.5  3.65  0.009
## 46  6 -833.383 1679.7  3.82  0.008
## 37  4 -835.626 1679.7  3.82  0.008
## 20  5 -834.554 1679.7  3.90  0.008
## 51  5 -834.573 1679.8  3.93  0.007
## 62  7 -832.416 1680.0  4.19  0.007
## 53  5 -834.715 1680.1  4.22  0.006
## 59  6 -833.597 1680.1  4.24  0.006
## 39  5 -834.746 1680.1  4.28  0.006
## 22  5 -834.747 1680.1  4.28  0.006
## 34  4 -835.855 1680.1  4.28  0.006
## 6   4 -835.901 1680.2  4.37  0.006
## 50  5 -834.834 1680.3  4.45  0.006
## 44  6 -833.875 1680.6  4.80  0.005
## 64  8 -831.590 1680.7  4.90  0.005
## 60  7 -832.841 1680.9  5.04  0.004
## 8   5 -835.167 1681.0  5.12  0.004
## 55  6 -834.153 1681.2  5.36  0.004
## 36  5 -835.366 1681.4  5.52  0.003
## 24  6 -834.284 1681.5  5.62  0.003
## 38  5 -835.595 1681.8  5.98  0.003
## 52  6 -834.550 1682.0  6.15  0.002
## 40  6 -834.645 1682.2  6.34  0.002
## 54  6 -834.713 1682.3  6.48  0.002
## 56  7 -834.140 1683.5  7.64  0.001
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | ID
out <- summary(model.avg(model_dredge, subset = delta <= 2))
out
## 
## Call:
## model.avg(object = model_dredge, subset = delta <= 2)
## 
## Component model call: 
## glmer(formula = round(RetreatHeight, 0) ~ <18 unique rhs>, data = 
##      webs_center, family = poisson, control = glmerControl(optimizer = 
##      "bobyqa", optCtrl = list(maxfun = 20000)))
## 
## Component models: 
##        df  logLik    AICc delta weight
## 4       3 -834.80 1675.84  0.00   0.09
## 5       3 -834.88 1676.00  0.16   0.08
## 45      4 -833.86 1676.14  0.30   0.08
## (Null)  2 -836.01 1676.15  0.30   0.08
## 234     5 -832.76 1676.15  0.31   0.08
## 34      4 -833.93 1676.27  0.43   0.07
## 345     5 -833.06 1676.75  0.91   0.06
## 14      4 -834.31 1677.04  1.20   0.05
## 24      4 -834.37 1677.15  1.31   0.05
## 2345    6 -832.18 1677.26  1.41   0.05
## 145     5 -833.33 1677.30  1.46   0.04
## 2       3 -835.53 1677.30  1.46   0.04
## 1234    6 -832.23 1677.35  1.50   0.04
## 134     5 -833.39 1677.42  1.58   0.04
## 25      4 -834.59 1677.60  1.76   0.04
## 2346    6 -832.38 1677.65  1.80   0.04
## 245     5 -833.60 1677.83  1.98   0.03
## 1345    6 -832.47 1677.84  2.00   0.03
## 
## Term codes: 
##                            log_patch                              log_tdh 
##                                    1                                    2 
##                              log_tdr                         log_TotalSub 
##                                    3                                    4 
## scale(road_length_m, center = FALSE)      scale(spec_rad, center = FALSE) 
##                                    5                                    6 
## 
## Model-averaged coefficients:  
## (full average) 
##                                      Estimate Std. Error Adjusted SE z value
## (Intercept)                           2.69800    2.74792     2.76264   0.977
## log_TotalSub                          0.58348    0.48668     0.48956   1.192
## scale(road_length_m, center = FALSE)  0.43258    0.70643     0.71051   0.609
## log_tdh                              -0.13792    0.26422     0.26586   0.519
## log_tdr                              -0.09461    0.14729     0.14799   0.639
## log_patch                             0.04317    0.12190     0.12273   0.352
## scale(spec_rad, center = FALSE)      -0.20758    1.60533     1.61706   0.128
##                                      Pr(>|z|)
## (Intercept)                             0.329
## log_TotalSub                            0.233
## scale(road_length_m, center = FALSE)    0.543
## log_tdh                                 0.604
## log_tdr                                 0.523
## log_patch                               0.725
## scale(spec_rad, center = FALSE)         0.898
##  
## (conditional average) 
##                                      Estimate Std. Error Adjusted SE z value
## (Intercept)                            2.6980     2.7479      2.7626   0.977
## log_TotalSub                           0.7725     0.4093      0.4138   1.867
## scale(road_length_m, center = FALSE)   1.0413     0.7533      0.7625   1.366
## log_tdh                               -0.3763     0.3175      0.3212   1.172
## log_tdr                               -0.2307     0.1466      0.1483   1.555
## log_patch                              0.2033     0.1934      0.1959   1.038
## scale(spec_rad, center = FALSE)       -5.6099     6.2722      6.3531   0.883
##                                      Pr(>|z|)  
## (Intercept)                            0.3288  
## log_TotalSub                           0.0619 .
## scale(road_length_m, center = FALSE)   0.1720  
## log_tdh                                0.2413  
## log_tdr                                0.1199  
## log_patch                              0.2994  
## scale(spec_rad, center = FALSE)        0.3772  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
info <- out[[10]]
#TotalSub
exp(info[2, 1])
## [1] 2.16528
exp(info[2, 2])
## [1] 1.505806
#road length
exp(info[3, 1])
## [1] 2.832813
exp(info[3, 2])
## [1] 2.123967
#tdh
exp(info[4, 1])
## [1] 0.6863857
exp(info[4, 2])
## [1] 1.373621
#tdr
exp(info[5, 1])
## [1] 0.7939575
exp(info[5, 2])
## [1] 1.157922
#patch
exp(info[6, 1])
## [1] 1.225422
exp(info[6, 2])
## [1] 1.213405
#spec_rad
exp(info[7, 1])
## [1] 0.00366157
exp(info[7, 2])
## [1] 529.6281
#car::vif(get.models(model_dredge, 1)[[1]]) # looking for less than 2

Overall Fixed effects: Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.1686 0.6215 5.098 3.43e-07 ** log_tdr -0.1896 0.1023 -1.854 0.06378 . negative* log_TotalSub 0.4812 0.1659 2.901 0.00372 ** positive

R-squared = 0.3671774 (marginal - fixed) 0.9521663 (conditional - fixed and random)

(conditional average) Estimate Std. Error Adjusted SE z value Pr(>|z|)
(Intercept) 3.3701 1.2954 1.3012 2.590 0.00959 ** log_tdr -0.1986 0.1022 0.1032 1.924 0.05432 .
log_TotalSub 0.4697 0.1814 0.1830 2.567 0.01026 * log_tdh -0.3307 0.2719 0.2745 1.205 0.22832

Urban Forest Fixed effects: Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.1283 0.8355 1.350 0.17687
tree100m_frac 5.0864 1.7266 2.946 0.00322 ** positive

R-squared = 0.4787743 (marginal - fixed) 0.9596909 (conditional - fixed and random)

(conditional average) Estimate Std. Error Adjusted SE z value Pr(>|z|)
(Intercept) 0.3764 1.3417 1.3881 0.271 0.7863
tree100m_frac 5.0033 1.7056 1.7823 2.807 0.0050 ** log_TotalSub 0.6815 0.3801 0.3976 1.714 0.0865 . log_tdr -0.2415 0.1683 0.1769 1.365 0.1721

Urban Center Fixed effects: Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.0705 0.4687 4.417 1e-05 ** log_TotalSub 0.6449 0.3921 1.645 0.1 positive*

R-squared = 0.1332201 (marginal - fixed) 0.9103755 (conditional - fixed and random)

(conditional average) Estimate Std. Error Adjusted SE z value Pr(>|z|) (Intercept) 2.6980 2.7479 2.7626 0.977 0.3288 log_TotalSub 0.7725 0.4093 0.4138 1.867 0.0619 . scale(road_length_m, center = FALSE) 1.0413 0.7533 0.7625 1.366 0.1720 log_tdh -0.3763 0.3175 0.3212 1.172 0.2413 log_tdr -0.2307 0.1466 0.1483 1.555 0.1199 log_patch 0.2033 0.1934 0.1959 1.038 0.2994 scale(spec_rad, center = FALSE) -5.6099 6.2722 6.3531 0.883 0.3772

colorc = "#d95f02" 
colorf = "#1b9e77"
coloro = "black"

# Plant Species Richness
predictions <- expand.grid(log_TotalSub = seq(0, 3.15, 0.03),
                           log_tdr = mean(webs$log_tdr))
predictions$response <- predict(top_model_wh_o, newdata = predictions, se.fit = TRUE, re.form = NA, type = "response")

myFunc <- function(mm) {
    predict(mm, newdata = predictions, re.form = ~0, type = "response")
}
#bigBoot_height_overall_plant <- bootMer(top_model_wh_o, myFunc, nsim = 1000)
#saveRDS(bigBoot_height_overall_plant, file = "bootstrapping/bigBoot_height_overall_plant.Rds")
bigBoot_height_overall_plant <- readRDS("bootstrapping/bigBoot_height_overall_plant.Rds")
predSE <- t(apply(bigBoot_height_overall_plant$t, MARGIN = 2, FUN = sd))
predictions$SE <- predSE[1, ]

height_plant <- ggplot(webs, aes(x = log_TotalSub, y = log(RetreatHeight))) + 
  geom_point(color = "grey", 
             size = 1) +
  geom_line(aes(x = log_TotalSub, y = log(response)), 
            size = 1, 
            color = coloro, 
            data = predictions) + 
  geom_ribbon(aes(x = log_TotalSub, 
                  ymin = log(response - 1.96 * SE), 
                  ymax = log(response + 1.96 * SE)), 
              fill = coloro, 
              data = predictions, 
              inherit.aes = FALSE, 
              alpha = 0.25) +
  scale_y_continuous(limits = c(min(log(webs$RetreatHeight)), max(log(webs$RetreatHeight)) + 0.1 * diff(range(log(webs$RetreatHeight))))) +  
  theme_classic() +
  xlab("Plant Species Richness [log()]") +
  ylab("Web Height [log(cm)]") + 
  theme(text = element_text(size = 10)) + 
  theme(axis.text.x = element_text(colour = "black", size = 10)) + 
  theme(axis.text.y = element_text(colour = "black", size = 10)) +
  theme(axis.text.x = element_blank(), 
        axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.y = element_blank()) +
  annotate(geom = "text", x = min(predictions$log_TotalSub) + diff(range(predictions$log_TotalSub))/2, y = max(log(webs$RetreatHeight)), label = "**", color = "black", size = 4)

# Road Disturbance
predictions <- expand.grid(log_TotalSub = mean(webs$log_TotalSub),
                           log_tdr = seq(2.73, 7.22, 0.05))
predictions$response <- predict(top_model_wh_o, newdata = predictions, se.fit = TRUE, re.form = NA, type = "response")

myFunc <- function(mm) {
    predict(mm, newdata = predictions, re.form = ~0, type = "response")
}
#bigBoot_height_overall_tdr <- bootMer(top_model_wh_o, myFunc, nsim = 1000)
#saveRDS(bigBoot_height_overall_tdr, file = "bootstrapping/bigBoot_height_overall_tdr.Rds")
bigBoot_height_overall_tdr <- readRDS("bootstrapping/bigBoot_height_overall_tdr.Rds")
predSE <- t(apply(bigBoot_height_overall_tdr$t, MARGIN = 2, FUN = sd))
predictions$SE <- predSE[1, ]

height_tdr <- ggplot(webs, aes(x = log_tdr, y = log(RetreatHeight))) + 
  geom_point(color = "grey", 
             size = 1) +
  geom_line(aes(x = log_tdr, y = log(response)), 
            size = 1, 
            color = coloro, 
            data = predictions) + 
  geom_ribbon(aes(x = log_tdr, 
                  ymin = log(response - 1.96 * SE), 
                  ymax = log(response + 1.96 * SE)), 
              fill = coloro, 
              data = predictions, 
              inherit.aes = FALSE, 
              alpha = 0.25) +
  scale_y_continuous(limits = c(min(log(webs$RetreatHeight)), max(log(webs$RetreatHeight)) + 0.1 * diff(range(log(webs$RetreatHeight))))) +  
  theme_classic() +
  xlab("Road Disturbance [log(vehicles/day/m)]") +
  ylab("Web Height [log(cm)]") + 
  theme(text = element_text(size = 10)) + 
  theme(axis.text.x = element_text(colour = "black", size = 10)) + 
  theme(axis.text.y = element_text(colour = "black", size = 10)) +
  theme(axis.text.x = element_blank(), 
        axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.y = element_blank()) +
  annotate(geom = "text", x = min(predictions$log_tdr) + diff(range(predictions$log_tdr))/2, y = max(log(webs$RetreatHeight)), label = "N. S.", color = "black", size = 4)

height_plant

height_tdr

colorc = "#d95f02" 
colorf = "#1b9e77"
coloro = "black"

# Tree
predictions <- expand.grid(tree100m_frac = seq(0.31, 0.66, 0.003))
predictions$response <- predict(top_model_wh_f, newdata = predictions, se.fit = TRUE, re.form = NA, type = "response")

myFunc <- function(mm) {
    predict(mm, newdata = predictions, re.form = ~0, type = "response")
}
#bigBoot_height_forest_tree <- bootMer(top_model_wh_f, myFunc, nsim = 1000)
#saveRDS(bigBoot_height_forest_tree, file = "bootstrapping/bigBoot_height_forest_tree.Rds")
bigBoot_height_forest_tree <- readRDS("bootstrapping/bigBoot_height_forest_tree.Rds")
predSE <- t(apply(bigBoot_height_forest_tree$t, MARGIN = 2, FUN = sd))
predictions$SE <- predSE[1, ]

height_forest_tree <- ggplot(webs_forest, aes(x = log(tree100m_frac), y = log(RetreatHeight))) + 
  geom_point(color = "grey", 
             size = 1) +
  geom_line(aes(x = log(tree100m_frac), y = log(response)), 
            size = 1, 
            color = colorf, 
            data = predictions) + 
  geom_ribbon(aes(x = log(tree100m_frac), 
                  ymin = log(response - 1.96 * SE), 
                  ymax = log(response + 1.96 * SE)), 
              fill = colorf, 
              data = predictions, 
              inherit.aes = FALSE, 
              alpha = 0.25) +
  scale_y_continuous(limits = c(min(log(webs_forest$RetreatHeight)), max(log(webs_forest$RetreatHeight)) + 0.1 * diff(range(log(webs_forest$RetreatHeight))))) +  
  theme_classic() +
  xlab("Proportion of Tree Cover in 100-Meter Buffer of Site [log()]") +
  ylab("Web Height [log(cm)]") + 
  theme(text = element_text(size = 10)) + 
  theme(axis.text.x = element_text(colour = "black", size = 10)) + 
  theme(axis.text.y = element_text(colour = "black", size = 10)) +
  theme(axis.text.x = element_blank(), 
        axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.y = element_blank()) +
  annotate(geom = "text", x = min(log(predictions$tree100m_frac)) + diff(range(log(predictions$tree100m_frac)))/2, y = max(log(webs_forest$RetreatHeight)), label = "**", color = "black", size = 4)

height_forest_tree

colorc = "#d95f02" 
colorf = "#1b9e77"
coloro = "black"

# Plant Species Richness
predictions <- expand.grid(log_TotalSub = seq(0, 1.80, 0.02))
predictions$response <- predict(top_model_wh_c, newdata = predictions, se.fit = TRUE, re.form = NA, type = "response")

myFunc <- function(mm) {
    predict(mm, newdata = predictions, re.form = ~0, type = "response")
}
#bigBoot_height_center_plant <- bootMer(top_model_wh_c, myFunc, nsim = 1000)
#saveRDS(bigBoot_height_center_plant, file = "bootstrapping/bigBoot_height_center_plant.Rds")
bigBoot_height_center_plant <- readRDS("bootstrapping/bigBoot_height_center_plant.Rds")
predSE <- t(apply(bigBoot_height_center_plant$t, MARGIN = 2, FUN = sd))
predictions$SE <- predSE[1, ]

height_center_plant <- ggplot(webs_center, aes(x = log_TotalSub, y = log(RetreatHeight))) + 
  geom_point(color = "grey", 
             size = 1) +
  geom_line(aes(x = log_TotalSub, y = log(response)), 
            size = 1, 
            color = colorc, 
            data = predictions) + 
  geom_ribbon(aes(x = log_TotalSub, 
                  ymin = log(response - 1.96 * SE), 
                  ymax = log(response + 1.96 * SE)), 
              fill = colorc, 
              data = predictions, 
              inherit.aes = FALSE, 
              alpha = 0.25) +
  scale_y_continuous(limits = c(min(log(webs_forest$RetreatHeight)) - 4, max(log(webs_forest$RetreatHeight)) + 0.1 * diff(range(log(webs_forest$RetreatHeight))))) +  
  theme_classic() +
  xlab("Plant Species Richness [log()]") +
  ylab("Web Height [log(cm)]") + 
  theme(text = element_text(size = 10)) + 
  theme(axis.text.x = element_text(colour = "black", size = 10)) + 
  theme(axis.text.y = element_text(colour = "black", size = 10)) +
  theme(axis.text.x = element_blank(), 
        axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        axis.title.y = element_blank()) +
  annotate(geom = "text", x = min((predictions$log_TotalSub)) + diff(range((predictions$log_TotalSub)))/2, y = max(log(webs_forest$RetreatHeight)), label = "N.S.", color = "black", size = 4)

height_center_plant

Web Height by Land

height_land <- glmer(round(RetreatHeight, 0) ~ Land + (1 | ID), data = webs, family = poisson)
summary(height_land)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: round(RetreatHeight, 0) ~ Land + (1 | ID)
##    Data: webs
## 
##      AIC      BIC   logLik deviance df.resid 
##   2457.7   2475.0  -1222.9   2445.7      125 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.4619 -2.2969 -0.5715  1.2929 23.5226 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 0.4474   0.6689  
## Number of obs: 131, groups:  ID, 22
## 
## Fixed effects:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)            2.7429     0.3387   8.100 5.52e-16 ***
## LandUrban, Medium     -0.2314     0.4395  -0.527   0.5984    
## LandUrban, Low         0.8390     0.5829   1.439   0.1501    
## LandDeciduous Forest   0.7705     0.4246   1.815   0.0696 .  
## LandWoody Wetlands     0.8142     0.5200   1.566   0.1174    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) LndU,M LndU,L LndDcF
## LndUrbn,Mdm -0.771                     
## LandUrbn,Lw -0.581  0.448              
## LndDcdsFrst -0.797  0.615  0.463       
## LndWdyWtlnd -0.651  0.502  0.378  0.519
car::Anova(height_land, type = 3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: round(RetreatHeight, 0)
##              Chisq Df Pr(>Chisq)    
## (Intercept) 65.603  1  5.516e-16 ***
## Land        10.495  4    0.03287 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
comp <- glht(height_land, linfct = mcp(Land = "Tukey"))
summary(comp)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: glmer(formula = round(RetreatHeight, 0) ~ Land + (1 | ID), data = webs, 
##     family = poisson)
## 
## Linear Hypotheses:
##                                        Estimate Std. Error z value Pr(>|z|)  
## Urban, Medium - Urban, High == 0       -0.23144    0.43946  -0.527   0.9842  
## Urban, Low - Urban, High == 0           0.83903    0.58292   1.439   0.5952  
## Deciduous Forest - Urban, High == 0     0.77051    0.42463   1.815   0.3581  
## Woody Wetlands - Urban, High == 0       0.81417    0.51995   1.566   0.5122  
## Urban, Low - Urban, Medium == 0         1.07047    0.55096   1.943   0.2882  
## Deciduous Forest - Urban, Medium == 0   1.00196    0.37955   2.640   0.0615 .
## Woody Wetlands - Urban, Medium == 0     1.04561    0.48379   2.161   0.1895  
## Deciduous Forest - Urban, Low == 0     -0.06852    0.53921  -0.127   0.9999  
## Woody Wetlands - Urban, Low == 0       -0.02486    0.61709  -0.040   1.0000  
## Woody Wetlands - Deciduous Forest == 0  0.04366    0.47039   0.093   1.0000  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
webs <- webs %>% 
  mutate(Land = fct_relevel(Land, "Deciduous Forest", 
                            "Woody Wetlands", 
                            "Urban, High", 
                            "Urban, Medium", 
                            "Urban, Low"))
predictions <- expand.grid(Land = levels(factor(webs$Land)))
predictions$response <- predict(height_land, newdata = predictions, se.fit = TRUE, re.form = NA, type = "response")

myFunc <- function(mm) {
    predict(mm, newdata = predictions, re.form = ~0, type = "response")
}
#bigBoot_height_land <- bootMer(height_land, myFunc, nsim = 1000)
#saveRDS(bigBoot_height_land, file = "bootstrapping/bigBoot_height_land.Rds")
bigBoot_height_land <- readRDS("bootstrapping/bigBoot_height_land.Rds")
predSE <- t(apply(bigBoot_height_land$t, MARGIN = 2, FUN = sd))
predictions$SE <- predSE[1, ]
#predictions2 <- summary(emmeans(height_land, ~Land,type = "response"))

ggplot(webs, aes(x = Land, y = RetreatHeight)) + 
  geom_jitter(color = "grey", width = 0.1, size = 0.5) +
  geom_point(aes(x = Land, y = response, color = Land), size = 1, data = predictions) + 
  geom_errorbar(aes(x = Land, 
                  ymin = response - SE, 
                  ymax = response + SE, color = Land), data = predictions, inherit.aes = FALSE, width = 0.25, size = 1) +
  theme_classic() +
  scale_color_manual(values = c("Deciduous Forest" = "#1b9e77", "Woody Wetlands" = "#1b9e77", "Urban, High" = "#d95f02", "Urban, Medium" = "#d95f02", "Urban, Low" = "#d95f02")) +
  xlab("Land Cover Class") +
  ylab("Web Height [log(cm)]") + 
  theme(text = element_text(size = 10)) + 
  theme(axis.text.x=element_text(colour="black", size=10)) + 
  theme(axis.text.y=element_text(colour="black", size=10)) +
  theme(legend.position = "none") +
  theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust=1)) +
  annotate(geom="text", x=1, y=max(webs$RetreatHeight), label="A", color="black", size = 3) +
  annotate(geom="text", x=2, y=max(webs$RetreatHeight), label="A", color="black", size = 3) +
  annotate(geom="text", x=3, y=max(webs$RetreatHeight), label="A", color="black", size = 3) +
  annotate(geom="text", x=4, y=max(webs$RetreatHeight), label="A", color="black", size = 3) +
  annotate(geom="text", x=5, y=max(webs$RetreatHeight), label="A", color="black", size = 3)

Web height varies significantly across land cover class (chi = 10.495, df = 4, p = 0.033).